Intro

Fit density (also referred to as CPUE) model with environmental predictors and use that to calculate weighted mean dissolved oxygen, temperature and depth of Baltic cod

# Load libraries, install if needed
library(tidyverse); theme_set(theme_light(base_size = 10))
library(readxl)
library(tidylog)
library(RCurl)
library(viridis)
library(RColorBrewer)
library(patchwork)
library(janitor)
library(icesDatras)
library(mapdata)
library(patchwork)
library(rgdal)
library(raster)
library(sf)
library(rgeos)
library(chron)
library(lattice)
library(ncdf4)
library(marmap)
library(rnaturalearth)
library(rnaturalearthdata)
library(mapplots)
library(EnvStats)
library(qwraps2)
#remotes::install_github("pbs-assess/sdmTMB")
library(sdmTMB)# To load entire cache in interactive r session, do: qwraps2::lazyload_cache_dir(path = "R/analysis/spatial_trend_models_cache/html")

For maps

# Specify map ranges
ymin = 55; ymax = 58; xmin = 14; xmax = 20

map_data <- rnaturalearth::ne_countries(
  scale = "medium",
  returnclass = "sf", continent = "europe")

# Crop the polygon for plotting and efficiency:
# st_bbox(map_data) # find the rough coordinates
swe_coast <- suppressWarnings(suppressMessages(
  st_crop(map_data,
          c(xmin = xmin, ymin = ymin, xmax = xmax, ymax = ymax))))

# Transform our map into UTM 33 coordinates, which is the equal-area projection we fit in:
utm_zone33 <- 32633
swe_coast_proj <- sf::st_transform(swe_coast, crs = utm_zone33)

ggplot(swe_coast_proj) + geom_sf() 


# Define plotting theme for main plot
theme_plot <- function(base_size = 10, base_family = "") {
  theme_light(base_size = 10, base_family = "") +
    theme(
      axis.text.x = element_text(angle = 90),
      axis.text = element_text(size = 8),
      legend.text = element_text(size = 8),
      legend.title = element_text(size = 8),
      legend.position = "bottom",
      legend.key.height = unit(0.2, "cm"),
      legend.margin = margin(0, 0, 0, 0),
      legend.box.margin = margin(-5, -5, -5, -5),
      strip.text = element_text(size = 8, colour = 'black', margin = margin()),
      strip.background = element_rect(fill = "grey90")
      )
}

# Define plotting theme for facet_wrap map with years
theme_facet_map <- function(base_size = 10, base_family = "") {
  theme_light(base_size = 10, base_family = "") +
    theme(
        axis.text.x = element_text(angle = 90),
        axis.text = element_text(size = 6),
        strip.text = element_text(size = 8, colour = 'black', margin = margin()),
        strip.background = element_rect(fill = "grey90")
      )
}

# Make default base map plot
plot_map_raster <- 
ggplot(swe_coast_proj) + 
  geom_sf(size = 0.3) +
  labs(x = "Longitude", y = "Latitude") +
  theme_facet_map(base_size = 14)

# Function to convert from lat lon to UTM
LongLatToUTM <- function(x, y, zone){
  xy <- data.frame(ID = 1:length(x), X = x, Y = y)
  coordinates(xy) <- c("X", "Y")
  proj4string(xy) <- CRS("+proj=longlat +datum=WGS84")  ## for example
  res <- spTransform(xy, CRS(paste("+proj=utm +zone=",zone," ellps=WGS84",sep='')))
  return(as.data.frame(res))
}

Read data

d <- readr::read_csv("/Users/maxlindmark/Desktop/R_STUDIO_PROJECTS/bentfish/data/for_analysis/stomach/full_stomach_data_21.10.25.csv") %>% dplyr::select(-X1)
#> Warning: Missing column names filled in: 'X1' [1]
#> 
#> ── Column specification ────────────────────────────────────────────────────────
#> cols(
#>   .default = col_double(),
#>   Sample = col_character(),
#>   Length.code = col_character(),
#>   Prey.sp.code = col_character(),
#>   Comment = col_character(),
#>   Parasites.in.stomach = col_character(),
#>   Coun. = col_character(),
#>   Ship = col_logical(),
#>   Cruise = col_character(),
#>   Predator.code = col_character(),
#>   Sex = col_character(),
#>   Q.year = col_character(),
#>   Date = col_date(format = ""),
#>   Index = col_character(),
#>   Ices.rect = col_character(),
#>   source = col_character(),
#>   Number = col_logical(),
#>   Sample.type = col_logical(),
#>   transect = col_logical(),
#>   Depthstep = col_logical(),
#>   Gonad.weight.roundfish = col_logical()
#>   # ... with 6 more columns
#> )
#> ℹ Use `spec()` for the full column specifications.
#> Warning: 94462 parsing failures.
#>  row  col           expected actual                                                                                                             file
#> 3364 Ship 1/0/T/F/TRUE/FALSE   SVEA '/Users/maxlindmark/Desktop/R_STUDIO_PROJECTS/bentfish/data/for_analysis/stomach/full_stomach_data_21.10.25.csv'
#> 3365 Ship 1/0/T/F/TRUE/FALSE   SVEA '/Users/maxlindmark/Desktop/R_STUDIO_PROJECTS/bentfish/data/for_analysis/stomach/full_stomach_data_21.10.25.csv'
#> 3366 Ship 1/0/T/F/TRUE/FALSE   SVEA '/Users/maxlindmark/Desktop/R_STUDIO_PROJECTS/bentfish/data/for_analysis/stomach/full_stomach_data_21.10.25.csv'
#> 3367 Ship 1/0/T/F/TRUE/FALSE   SVEA '/Users/maxlindmark/Desktop/R_STUDIO_PROJECTS/bentfish/data/for_analysis/stomach/full_stomach_data_21.10.25.csv'
#> 3368 Ship 1/0/T/F/TRUE/FALSE   SVEA '/Users/maxlindmark/Desktop/R_STUDIO_PROJECTS/bentfish/data/for_analysis/stomach/full_stomach_data_21.10.25.csv'
#> .... .... .................. ...... ................................................................................................................
#> See problems(...) for more details.

# Add UTM coordinates
utm_coords <- LongLatToUTM(d$Long, d$Lat, zone = 33)
#> Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
#> prefer_proj): Discarded datum Unknown based on WGS84 ellipsoid in CRS definition
d$X <- utm_coords$X/1000
d$Y <- utm_coords$Y/1000

Explore data

head(data.frame(d))
#>   SubDiv Year Month Day N.with.food N.regurgitated N.skeletal N.empty HN Sample
#> 1     28 2017    11  27           1             NA         NA      NA 52    507
#> 2     28 2017    11  27           1             NA         NA      NA 52    507
#> 3     28 2017    11  27           1             NA         NA      NA 52    507
#> 4     28 2017    11  27           1             NA         NA      NA 52    507
#> 5     28 2017    11  27           1             NA         NA      NA 52    507
#> 6     28 2017    11  27           1             NA         NA      NA 52    507
#>   Length.code    Prey.sp.code Prey.size Prey.weight Prey.nr Stage.digestion
#> 1        <NA> Saduria entomon        20        0.21       1               1
#> 2        <NA> Saduria entomon        16        0.10       1               1
#> 3        <NA> Saduria entomon        20        0.15       1               1
#> 4        <NA> Saduria entomon        20        0.17       1               1
#> 5        <NA> Saduria entomon        NA        0.27       2               1
#> 6        <NA> Macoma balthica        NA        0.19       1               1
#>   Comment Parasites.in.stomach Quarter Coun. Ship Cruise Pred.size.mm
#> 1    <NA>                 <NA>       4   SWE   NA   BITS          200
#> 2    <NA>                 <NA>       4   SWE   NA   BITS          200
#> 3    <NA>                 <NA>       4   SWE   NA   BITS          200
#> 4    <NA>                 <NA>       4   SWE   NA   BITS          200
#> 5    <NA>                 <NA>       4   SWE   NA   BITS          200
#> 6    <NA>                 <NA>       4   SWE   NA   BITS          200
#>   Pred.weight.g Predator.gutted.weight Predator.code Sex Maturity Gall.content
#> 1           100                     92           FLE   M        4           NA
#> 2           100                     92           FLE   M        4           NA
#> 3           100                     92           FLE   M        4           NA
#> 4           100                     92           FLE   M        4           NA
#> 5           100                     92           FLE   M        4           NA
#> 6           100                     92           FLE   M        4           NA
#>   Q.year Age       Date        Index      Lat     Long Depth.catch Ices.rect
#> 1 4_2017   3 2017-11-27 OXBH.2017.52 57.46667 19.23333          67      43G9
#> 2 4_2017   3 2017-11-27 OXBH.2017.52 57.46667 19.23333          67      43G9
#> 3 4_2017   3 2017-11-27 OXBH.2017.52 57.46667 19.23333          67      43G9
#> 4 4_2017   3 2017-11-27 OXBH.2017.52 57.46667 19.23333          67      43G9
#> 5 4_2017   3 2017-11-27 OXBH.2017.52 57.46667 19.23333          67      43G9
#> 6 4_2017   3 2017-11-27 OXBH.2017.52 57.46667 19.23333          67      43G9
#>        source Number Sample.type transect Depthstep Gonad.weight.roundfish
#> 1 Max Postdoc     NA          NA       NA        NA                     NA
#> 2 Max Postdoc     NA          NA       NA        NA                     NA
#> 3 Max Postdoc     NA          NA       NA        NA                     NA
#> 4 Max Postdoc     NA          NA       NA        NA                     NA
#> 5 Max Postdoc     NA          NA       NA        NA                     NA
#> 6 Max Postdoc     NA          NA       NA        NA                     NA
#>   Liver.weight.roundfish Stomach.content Perc.stomac.content comments Processed
#> 1                     NA              NA                  NA       NA        NA
#> 2                     NA              NA                  NA       NA        NA
#> 3                     NA              NA                  NA       NA        NA
#> 4                     NA              NA                  NA       NA        NA
#> 5                     NA              NA                  NA       NA        NA
#> 6                     NA              NA                  NA       NA        NA
#>      Unique.pred.id       X        Y
#> 1 2017_4_52_507_FLE 753.841 6377.247
#> 2 2017_4_52_507_FLE 753.841 6377.247
#> 3 2017_4_52_507_FLE 753.841 6377.247
#> 4 2017_4_52_507_FLE 753.841 6377.247
#> 5 2017_4_52_507_FLE 753.841 6377.247
#> 6 2017_4_52_507_FLE 753.841 6377.247
sort(colnames(d))
#>  [1] "Age"                    "Comment"                "comments"              
#>  [4] "Coun."                  "Cruise"                 "Date"                  
#>  [7] "Day"                    "Depth.catch"            "Depthstep"             
#> [10] "Gall.content"           "Gonad.weight.roundfish" "HN"                    
#> [13] "Ices.rect"              "Index"                  "Lat"                   
#> [16] "Length.code"            "Liver.weight.roundfish" "Long"                  
#> [19] "Maturity"               "Month"                  "N.empty"               
#> [22] "N.regurgitated"         "N.skeletal"             "N.with.food"           
#> [25] "Number"                 "Parasites.in.stomach"   "Perc.stomac.content"   
#> [28] "Pred.size.mm"           "Pred.weight.g"          "Predator.code"         
#> [31] "Predator.gutted.weight" "Prey.nr"                "Prey.size"             
#> [34] "Prey.sp.code"           "Prey.weight"            "Processed"             
#> [37] "Q.year"                 "Quarter"                "Sample"                
#> [40] "Sample.type"            "Sex"                    "Ship"                  
#> [43] "source"                 "Stage.digestion"        "Stomach.content"       
#> [46] "SubDiv"                 "transect"               "Unique.pred.id"        
#> [49] "X"                      "Y"                      "Year"
#  [1] "Age"                    "Comment"                "comments"               "Coun."                 
#  [5] "Cruise"                 "Date"                   "Day"                    "Depth.catch"           
#  [9] "Depthstep"              "Gall.content"           "Gonad.weight.roundfish" "HN"                    
# [13] "Ices.rect"              "Index"                  "Lat"                    "Length.code"           
# [17] "Liver.weight.roundfish" "Long"                   "Maturity"               "Month"                 
# [21] "N.empty"                "N.regurgitated"         "N.skeletal"             "N.with.food"           
# [25] "Number"                 "Parasites.in.stomach"   "Perc.stomac.content"    "Pred.size.mm"          
# [29] "Pred.weight.g"          "Predator.code"          "Predator.gutted.weight" "Prey.nr"               
# [33] "Prey.size"              "Prey.sp.code"           "Prey.weight"            "Processed"             
# [37] "Q.year"                 "Quarter"                "Sample"                 "Sample.type"           
# [41] "Sex"                    "Ship"                   "source"                 "Stage.digestion"       
# [45] "Stomach.content"        "SubDiv"                 "transect"               "Unique.pred.id"        
# [49] "X"

# Plot data in space, color by survey
plot_map_raster +
  geom_point(data = d, aes(x = X * 1000, y = Y * 1000, color = Cruise), size = 0.5) +
  scale_color_brewer(palette = "Set2") + 
  facet_wrap(~ Cruise, ncol = 2) + 
  theme_plot()


# Calculate total weight of prey by predator ID and prey species (i.e., across prey sizes)
d_wide <- d %>% 
  drop_na(Prey.weight) %>% 
  group_by(Unique.pred.id, Prey.sp.code) %>% 
  summarise(tot_biom_per_prey = sum(Prey.weight)) %>% 
  ungroup() %>% 
  pivot_wider(names_from = Prey.sp.code, values_from = tot_biom_per_prey) %>% 
  mutate(across(everything(), ~replace_na(.x, 0))) %>% # Replace NA with 0
  janitor::clean_names()
#> drop_na: removed 1,216 rows (6%), 20,515 rows remaining
#> group_by: 2 grouping variables (Unique.pred.id, Prey.sp.code)
#> summarise: now 16,767 rows and 3 columns, one group variable remaining (Unique.pred.id)
#> ungroup: no grouping variables
#> pivot_wider: reorganized (Prey.sp.code, tot_biom_per_prey) into (Sprattus sprattus, Bylgides sarsi, Saduria entomon, Clupea harengus, Pisces, …) [was 16767x3, now 8311x113]
#> mutate: changed 7,085 values (85%) of 'Sprattus sprattus' (7085 fewer NA)
#>         changed 7,030 values (85%) of 'Bylgides sarsi' (7030 fewer NA)
#>         changed 6,431 values (77%) of 'Saduria entomon' (6431 fewer NA)
#>         changed 7,616 values (92%) of 'Clupea harengus' (7616 fewer NA)
#>         changed 7,380 values (89%) of 'Pisces' (7380 fewer NA)
#>         changed 8,185 values (98%) of 'Gadus morhua' (8185 fewer NA)
#>         changed 7,956 values (96%) of 'Cumacea' (7956 fewer NA)
#>         changed 7,345 values (88%) of 'Limecola balthica' (7345 fewer NA)
#>         changed 6,845 values (82%) of 'Mysis mixta' (6845 fewer NA)
#>         changed 8,307 values (>99%) of 'Corophium volutator' (8307 fewer NA)
#>         changed 7,639 values (92%) of 'Halicryptus spinulosus' (7639 fewer NA)
#>         changed 8,203 values (99%) of 'Stone' (8203 fewer NA)
#>         changed 7,662 values (92%) of 'Pontoporeia femorata' (7662 fewer NA)
#>         changed 8,285 values (>99%) of 'Enchelyopus cimbrius' (8285 fewer NA)
#>         changed 7,946 values (96%) of 'Clupeidae' (7946 fewer NA)
#>         changed 7,916 values (95%) of 'Gobiidae' (7916 fewer NA)
#>         changed 8,109 values (98%) of 'Neomysis integer' (8109 fewer NA)
#>         changed 7,989 values (96%) of 'Monoporeia affinis' (7989 fewer NA)
#>         changed 7,865 values (95%) of 'Gammarus sp.' (7865 fewer NA)
#>         changed 8,034 values (97%) of 'Mysidae' (8034 fewer NA)
#>         changed 8,265 values (99%) of 'Scoloplos armiger' (8265 fewer NA)
#>         changed 8,297 values (>99%) of 'plastic' (8297 fewer NA)
#>         changed 8,216 values (99%) of 'Priapulus caudatus' (8216 fewer NA)
#>         changed 8,257 values (99%) of 'Bivalvia' (8257 fewer NA)
#>         changed 8,295 values (>99%) of 'Hyperia galba' (8295 fewer NA)
#>         changed 8,091 values (97%) of 'Mytilus sp.' (8091 fewer NA)
#>         changed 8,019 values (96%) of 'Crangon crangon' (8019 fewer NA)
#>         changed 8,263 values (99%) of 'Idotea balthica' (8263 fewer NA)
#>         changed 8,309 values (>99%) of 'Trachurus trachurus' (8309 fewer NA)
#>         changed 8,063 values (97%) of 'Gasterosteus aculeatus' (8063 fewer NA)
#>         changed 8,309 values (>99%) of 'Annelida' (8309 fewer NA)
#>         changed 8,221 values (99%) of 'Algae' (8221 fewer NA)
#>         changed 8,296 values (>99%) of 'Platichthys flesus' (8296 fewer NA)
#>         changed 8,307 values (>99%) of 'Priapulidae' (8307 fewer NA)
#>         changed 8,306 values (>99%) of 'Waste' (8306 fewer NA)
#>         changed 8,278 values (>99%) of 'Praunus flexuosus' (8278 fewer NA)
#>         changed 8,301 values (>99%) of 'Unidentified mass' (8301 fewer NA)
#>         changed 8,308 values (>99%) of 'Merlangius merlangus' (8308 fewer NA)
#>         changed 8,252 values (99%) of 'Scales' (8252 fewer NA)
#>         changed 8,308 values (>99%) of 'Gadidae' (8308 fewer NA)
#>         changed 8,257 values (99%) of 'Crustacea' (8257 fewer NA)
#>         changed 8,211 values (99%) of 'Amphipoda' (8211 fewer NA)
#>         changed 8,258 values (99%) of 'Spine' (8258 fewer NA)
#>         changed 8,309 values (>99%) of 'Pleuronectes platessa' (8309 fewer NA)
#>         changed 8,310 values (>99%) of 'Anguilla anguilla' (8310 fewer NA)
#>         changed 8,310 values (>99%) of 'Empty' (8310 fewer NA)
#>         changed 8,310 values (>99%) of 'Mollusca' (8310 fewer NA)
#>         changed 8,309 values (>99%) of 'Pungitius pungitius' (8309 fewer NA)
#>         changed 7,462 values (90%) of 'Diastylis rathkei' (7462 fewer NA)
#>         changed 8,310 values (>99%) of 'Terebellides stroemii' (8310 fewer NA)
#>         changed 8,293 values (>99%) of 'Palaemon sp.' (8293 fewer NA)
#>         changed 8,305 values (>99%) of 'Ammodytes tobianus' (8305 fewer NA)
#>         changed 7,506 values (90%) of 'NA' (7506 fewer NA)
#>         changed 8,309 values (>99%) of 'Cottidae' (8309 fewer NA)
#>         changed 8,303 values (>99%) of 'Hediste diversicolor' (8303 fewer NA)
#>         changed 8,296 values (>99%) of 'Palaemon elegans' (8296 fewer NA)
#>         changed 8,308 values (>99%) of 'Spinachia spinachia' (8308 fewer NA)
#>         changed 8,304 values (>99%) of 'Zoarces viviparus' (8304 fewer NA)
#>         changed 8,286 values (>99%) of 'Ammodytidae' (8286 fewer NA)
#>         changed 8,304 values (>99%) of 'Caridea' (8304 fewer NA)
#>         changed 8,310 values (>99%) of 'Amphibalanus improvisus' (8310 fewer NA)
#>         changed 8,310 values (>99%) of 'Myoxocephalus quadricornis' (8310 fewer NA)
#>         changed 8,305 values (>99%) of 'Phyllodocida' (8305 fewer NA)
#>         changed 8,226 values (99%) of 'Remains' (8226 fewer NA)
#>         changed 8,302 values (>99%) of 'Palaemonidae' (8302 fewer NA)
#>         changed 8,298 values (>99%) of 'Carcinus maenas' (8298 fewer NA)
#>         changed 8,310 values (>99%) of 'Gastropoda' (8310 fewer NA)
#>         changed 8,300 values (>99%) of 'Sand' (8300 fewer NA)
#>         changed 8,308 values (>99%) of 'Cerastoderma glaucum' (8308 fewer NA)
#>         changed 8,302 values (>99%) of 'Hyperoplus lanceolatus' (8302 fewer NA)
#>         changed 8,285 values (>99%) of 'Mya arenaria' (8285 fewer NA)
#>         changed 8,298 values (>99%) of 'Hydrobia sp.' (8298 fewer NA)
#>         changed 8,304 values (>99%) of 'Wood' (8304 fewer NA)
#>         changed 8,308 values (>99%) of 'Pleuronectiformes' (8308 fewer NA)
#>         changed 8,310 values (>99%) of 'Scophthalmus maximus' (8310 fewer NA)
#>         changed 8,283 values (>99%) of 'Polychaeta' (8283 fewer NA)
#>         changed 8,213 values (99%) of 'Priapulida' (8213 fewer NA)
#>         changed 8,310 values (>99%) of 'Halicryptus' (8310 fewer NA)
#>         changed 8,309 values (>99%) of 'Neogobius melanostomus' (8309 fewer NA)
#>         changed 8,309 values (>99%) of 'Gobius niger' (8309 fewer NA)
#>         changed 8,310 values (>99%) of 'digestive tract' (8310 fewer NA)
#>         changed 8,310 values (>99%) of 'Pleuronectidae' (8310 fewer NA)
#>         changed 8,310 values (>99%) of 'sprattus sprattus' (8310 fewer NA)
#>         changed 8,310 values (>99%) of 'Gasterosteidae' (8310 fewer NA)
#>         changed 8,308 values (>99%) of 'litter' (8308 fewer NA)
#>         changed 8,289 values (>99%) of 'Mysida' (8289 fewer NA)
#>         changed 8,310 values (>99%) of 'Carbon' (8310 fewer NA)
#>         changed 8,282 values (>99%) of 'Copepoda' (8282 fewer NA)
#>         changed 8,309 values (>99%) of 'Calanoida' (8309 fewer NA)
#>         changed 8,310 values (>99%) of 'Belone belone' (8310 fewer NA)
#>         changed 8,212 values (99%) of 'Pontoporeiidae' (8212 fewer NA)
#>         changed 8,307 values (>99%) of 'Mucus' (8307 fewer NA)
#>         changed 8,310 values (>99%) of 'Pectinaria sp.' (8310 fewer NA)
#>         changed 8,191 values (99%) of 'Macoma balthica' (8191 fewer NA)
#>         changed 8,026 values (97%) of 'remains' (8026 fewer NA)
#>         changed 8,258 values (99%) of 'stone' (8258 fewer NA)
#>         changed 8,310 values (>99%) of 'carbon' (8310 fewer NA)
#>         changed 8,309 values (>99%) of 'Nephtys ciliata' (8309 fewer NA)
#>         changed 8,310 values (>99%) of 'Gastrosacus' (8310 fewer NA)
#>         changed 8,310 values (>99%) of 'wood' (8310 fewer NA)
#>         changed 8,310 values (>99%) of 'Litter/plastics' (8310 fewer NA)
#>         changed 8,310 values (>99%) of 'clupeidae' (8310 fewer NA)
#>         changed 8,298 values (>99%) of 'Pontoporeidae' (8298 fewer NA)
#>         changed 8,304 values (>99%) of 'sand' (8304 fewer NA)
#>         changed 8,309 values (>99%) of 'Decapoda' (8309 fewer NA)
#>         changed 8,310 values (>99%) of 'Agonus cataphractus' (8310 fewer NA)
#>         changed 8,310 values (>99%) of 'clupea harengus' (8310 fewer NA)
#>         changed 8,295 values (>99%) of 'Halicryptus spinolusus' (8295 fewer NA)
#>         changed 8,156 values (98%) of 'Mytilus edulis' (8156 fewer NA)
#>         changed 8,303 values (>99%) of 'priapulida' (8303 fewer NA)
#>         changed 8,310 values (>99%) of 'Prapulida' (8310 fewer NA)
#>         changed 8,310 values (>99%) of 'Myoxocephalus scorpius' (8310 fewer NA)

head(d_wide)
#> # A tibble: 6 x 113
#>   unique_pred_id sprattus_spratt… bylgides_sarsi saduria_entomon clupea_harengus
#>   <chr>                     <dbl>          <dbl>           <dbl>           <dbl>
#> 1 2007_4_101_20…             5.49           0                0               0  
#> 2 2007_4_101_20…             0              0.24             0               0  
#> 3 2007_4_101_20…             0              0                0.9             0  
#> 4 2007_4_101_20…            26.6            0                0             180. 
#> 5 2007_4_101_20…             0              0                0              88.5
#> 6 2007_4_101_20…            27.3            0                0             163. 
#> # … with 108 more variables: pisces <dbl>, gadus_morhua <dbl>, cumacea <dbl>,
#> #   limecola_balthica <dbl>, mysis_mixta <dbl>, corophium_volutator <dbl>,
#> #   halicryptus_spinulosus <dbl>, stone <dbl>, pontoporeia_femorata <dbl>,
#> #   enchelyopus_cimbrius <dbl>, clupeidae <dbl>, gobiidae <dbl>,
#> #   neomysis_integer <dbl>, monoporeia_affinis <dbl>, gammarus_sp <dbl>,
#> #   mysidae <dbl>, scoloplos_armiger <dbl>, plastic <dbl>,
#> #   priapulus_caudatus <dbl>, bivalvia <dbl>, hyperia_galba <dbl>,
#> #   mytilus_sp <dbl>, crangon_crangon <dbl>, idotea_balthica <dbl>,
#> #   trachurus_trachurus <dbl>, gasterosteus_aculeatus <dbl>, annelida <dbl>,
#> #   algae <dbl>, platichthys_flesus <dbl>, priapulidae <dbl>, waste <dbl>,
#> #   praunus_flexuosus <dbl>, unidentified_mass <dbl>,
#> #   merlangius_merlangus <dbl>, scales <dbl>, gadidae <dbl>, crustacea <dbl>,
#> #   amphipoda <dbl>, spine <dbl>, pleuronectes_platessa <dbl>,
#> #   anguilla_anguilla <dbl>, empty <dbl>, mollusca <dbl>,
#> #   pungitius_pungitius <dbl>, diastylis_rathkei <dbl>,
#> #   terebellides_stroemii <dbl>, palaemon_sp <dbl>, ammodytes_tobianus <dbl>,
#> #   na <dbl>, cottidae <dbl>, hediste_diversicolor <dbl>,
#> #   palaemon_elegans <dbl>, spinachia_spinachia <dbl>, zoarces_viviparus <dbl>,
#> #   ammodytidae <dbl>, caridea <dbl>, amphibalanus_improvisus <dbl>,
#> #   myoxocephalus_quadricornis <dbl>, phyllodocida <dbl>, remains <dbl>,
#> #   palaemonidae <dbl>, carcinus_maenas <dbl>, gastropoda <dbl>, sand <dbl>,
#> #   cerastoderma_glaucum <dbl>, hyperoplus_lanceolatus <dbl>,
#> #   mya_arenaria <dbl>, hydrobia_sp <dbl>, wood <dbl>, pleuronectiformes <dbl>,
#> #   scophthalmus_maximus <dbl>, polychaeta <dbl>, priapulida <dbl>,
#> #   halicryptus <dbl>, neogobius_melanostomus <dbl>, gobius_niger <dbl>,
#> #   digestive_tract <dbl>, pleuronectidae <dbl>, sprattus_sprattus_2 <dbl>,
#> #   gasterosteidae <dbl>, litter <dbl>, mysida <dbl>, carbon <dbl>,
#> #   copepoda <dbl>, calanoida <dbl>, belone_belone <dbl>, pontoporeiidae <dbl>,
#> #   mucus <dbl>, pectinaria_sp <dbl>, macoma_balthica <dbl>, remains_2 <dbl>,
#> #   stone_2 <dbl>, carbon_2 <dbl>, nephtys_ciliata <dbl>, gastrosacus <dbl>,
#> #   wood_2 <dbl>, litter_plastics <dbl>, clupeidae_2 <dbl>,
#> #   pontoporeidae <dbl>, sand_2 <dbl>, …
str(d_wide)
#> tibble [8,311 × 113] (S3: tbl_df/tbl/data.frame)
#>  $ unique_pred_id            : chr [1:8311] "2007_4_101_2007_11_7_101_DAN2_40G5_1_M_490_1028_1_COD" "2007_4_101_2007_11_7_101_DAN2_40G5_118_F_210_82_1_COD" "2007_4_101_2007_11_7_101_DAN2_40G5_45_F_480_784_1_COD" "2007_4_101_2007_11_7_101_DAN2_40G5_65_F_920_7140_1_COD" ...
#>  $ sprattus_sprattus         : num [1:8311] 5.49 0 0 26.59 0 ...
#>  $ bylgides_sarsi            : num [1:8311] 0 0.24 0 0 0 0 0 0 0 0 ...
#>  $ saduria_entomon           : num [1:8311] 0 0 0.9 0 0 0 0 0 0 0 ...
#>  $ clupea_harengus           : num [1:8311] 0 0 0 179.7 88.5 ...
#>  $ pisces                    : num [1:8311] 0 0 0 20.8 0 ...
#>  $ gadus_morhua              : num [1:8311] 0 0 0 0 114 ...
#>  $ cumacea                   : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ limecola_balthica         : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ mysis_mixta               : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ corophium_volutator       : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ halicryptus_spinulosus    : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ stone                     : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ pontoporeia_femorata      : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ enchelyopus_cimbrius      : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ clupeidae                 : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ gobiidae                  : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ neomysis_integer          : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ monoporeia_affinis        : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ gammarus_sp               : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ mysidae                   : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ scoloplos_armiger         : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ plastic                   : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ priapulus_caudatus        : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ bivalvia                  : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ hyperia_galba             : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ mytilus_sp                : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ crangon_crangon           : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ idotea_balthica           : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ trachurus_trachurus       : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ gasterosteus_aculeatus    : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ annelida                  : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ algae                     : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ platichthys_flesus        : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ priapulidae               : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ waste                     : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ praunus_flexuosus         : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ unidentified_mass         : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ merlangius_merlangus      : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ scales                    : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ gadidae                   : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ crustacea                 : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ amphipoda                 : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ spine                     : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ pleuronectes_platessa     : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ anguilla_anguilla         : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ empty                     : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ mollusca                  : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ pungitius_pungitius       : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ diastylis_rathkei         : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ terebellides_stroemii     : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ palaemon_sp               : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ ammodytes_tobianus        : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ na                        : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ cottidae                  : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ hediste_diversicolor      : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ palaemon_elegans          : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ spinachia_spinachia       : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ zoarces_viviparus         : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ ammodytidae               : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ caridea                   : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ amphibalanus_improvisus   : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ myoxocephalus_quadricornis: num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ phyllodocida              : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ remains                   : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ palaemonidae              : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ carcinus_maenas           : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ gastropoda                : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ sand                      : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ cerastoderma_glaucum      : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ hyperoplus_lanceolatus    : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ mya_arenaria              : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ hydrobia_sp               : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ wood                      : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ pleuronectiformes         : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ scophthalmus_maximus      : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ polychaeta                : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ priapulida                : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ halicryptus               : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ neogobius_melanostomus    : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ gobius_niger              : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ digestive_tract           : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ pleuronectidae            : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ sprattus_sprattus_2       : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ gasterosteidae            : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ litter                    : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ mysida                    : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ carbon                    : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ copepoda                  : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ calanoida                 : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ belone_belone             : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ pontoporeiidae            : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ mucus                     : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ pectinaria_sp             : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ macoma_balthica           : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ remains_2                 : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ stone_2                   : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ carbon_2                  : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ nephtys_ciliata           : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>   [list output truncated]
colnames(d_wide)
#>   [1] "unique_pred_id"             "sprattus_sprattus"         
#>   [3] "bylgides_sarsi"             "saduria_entomon"           
#>   [5] "clupea_harengus"            "pisces"                    
#>   [7] "gadus_morhua"               "cumacea"                   
#>   [9] "limecola_balthica"          "mysis_mixta"               
#>  [11] "corophium_volutator"        "halicryptus_spinulosus"    
#>  [13] "stone"                      "pontoporeia_femorata"      
#>  [15] "enchelyopus_cimbrius"       "clupeidae"                 
#>  [17] "gobiidae"                   "neomysis_integer"          
#>  [19] "monoporeia_affinis"         "gammarus_sp"               
#>  [21] "mysidae"                    "scoloplos_armiger"         
#>  [23] "plastic"                    "priapulus_caudatus"        
#>  [25] "bivalvia"                   "hyperia_galba"             
#>  [27] "mytilus_sp"                 "crangon_crangon"           
#>  [29] "idotea_balthica"            "trachurus_trachurus"       
#>  [31] "gasterosteus_aculeatus"     "annelida"                  
#>  [33] "algae"                      "platichthys_flesus"        
#>  [35] "priapulidae"                "waste"                     
#>  [37] "praunus_flexuosus"          "unidentified_mass"         
#>  [39] "merlangius_merlangus"       "scales"                    
#>  [41] "gadidae"                    "crustacea"                 
#>  [43] "amphipoda"                  "spine"                     
#>  [45] "pleuronectes_platessa"      "anguilla_anguilla"         
#>  [47] "empty"                      "mollusca"                  
#>  [49] "pungitius_pungitius"        "diastylis_rathkei"         
#>  [51] "terebellides_stroemii"      "palaemon_sp"               
#>  [53] "ammodytes_tobianus"         "na"                        
#>  [55] "cottidae"                   "hediste_diversicolor"      
#>  [57] "palaemon_elegans"           "spinachia_spinachia"       
#>  [59] "zoarces_viviparus"          "ammodytidae"               
#>  [61] "caridea"                    "amphibalanus_improvisus"   
#>  [63] "myoxocephalus_quadricornis" "phyllodocida"              
#>  [65] "remains"                    "palaemonidae"              
#>  [67] "carcinus_maenas"            "gastropoda"                
#>  [69] "sand"                       "cerastoderma_glaucum"      
#>  [71] "hyperoplus_lanceolatus"     "mya_arenaria"              
#>  [73] "hydrobia_sp"                "wood"                      
#>  [75] "pleuronectiformes"          "scophthalmus_maximus"      
#>  [77] "polychaeta"                 "priapulida"                
#>  [79] "halicryptus"                "neogobius_melanostomus"    
#>  [81] "gobius_niger"               "digestive_tract"           
#>  [83] "pleuronectidae"             "sprattus_sprattus_2"       
#>  [85] "gasterosteidae"             "litter"                    
#>  [87] "mysida"                     "carbon"                    
#>  [89] "copepoda"                   "calanoida"                 
#>  [91] "belone_belone"              "pontoporeiidae"            
#>  [93] "mucus"                      "pectinaria_sp"             
#>  [95] "macoma_balthica"            "remains_2"                 
#>  [97] "stone_2"                    "carbon_2"                  
#>  [99] "nephtys_ciliata"            "gastrosacus"               
#> [101] "wood_2"                     "litter_plastics"           
#> [103] "clupeidae_2"                "pontoporeidae"             
#> [105] "sand_2"                     "decapoda"                  
#> [107] "agonus_cataphractus"        "clupea_harengus_2"         
#> [109] "halicryptus_spinolusus"     "mytilus_edulis"            
#> [111] "priapulida_2"               "prapulida"                 
#> [113] "myoxocephalus_scorpius"
 
# Now make some calculations and aggregate some taxonomic level
d_wide2 <- d_wide %>% 
  mutate(amphipoda_tot = bylgides_sarsi + hyperia_galba + gammarus_sp + monoporeia_affinis + 
           corophium_volutator + amphipoda,
         bivalvia_tot = bivalvia + mytilus_sp + cerastoderma_glaucum + mya_arenaria + macoma_balthica + 
           mytilus_edulis,
         clupeidae_tot = clupeidae + clupea_harengus_2 + sprattus_sprattus_2 + clupeidae_2,
         sprattus_sprattus_tot = sprattus_sprattus + sprattus_sprattus_2,
         clupea_harengus_tot = clupea_harengus + clupea_harengus_2,
         #gadus_morhua_tot = gadus_morhua,
         gadiformes_tot = gadidae + merlangius_merlangus,
         gobiidae_tot = gobiidae,
         mysidae_tot = mysidae + neomysis_integer + mysis_mixta + mysida + gastrosacus,
         non_bio_tot = stone + plastic + sand + wood + litter + carbon + stone_2 + carbon_2 + wood_2 + 
           litter_plastics + sand_2,
         other_crustacea_tot = pontoporeia_femorata + crangon_crangon + idotea_balthica + 
           praunus_flexuosus + crustacea + diastylis_rathkei + palaemon_sp + palaemon_elegans + caridea +
           amphibalanus_improvisus + palaemonidae + carcinus_maenas + copepoda + calanoida + pontoporeiidae + decapoda, 
         other_tot = halicryptus_spinulosus + priapulus_caudatus + annelida + algae + priapulidae + waste +
           unidentified_mass + spine + empty + mollusca + na + remains + gastropoda + hydrobia_sp + 
           priapulida + halicryptus + digestive_tract + mucus + remains_2 + pontoporeidae +
           halicryptus_spinolusus + priapulida_2 + prapulida,
         other_pisces_tot = enchelyopus_cimbrius + trachurus_trachurus + gasterosteus_aculeatus + 
           scales + pleuronectes_platessa + anguilla_anguilla + pungitius_pungitius + ammodytes_tobianus +
           cottidae + spinachia_spinachia + zoarces_viviparus + ammodytidae + myoxocephalus_quadricornis +
           hyperoplus_lanceolatus + pleuronectiformes + scophthalmus_maximus + neogobius_melanostomus + gobius_niger +
           pleuronectidae + gasterosteidae + belone_belone + agonus_cataphractus + myoxocephalus_scorpius,
         platichthys_flesus_tot = platichthys_flesus,
         polychaeta_tot = bylgides_sarsi + scoloplos_armiger + terebellides_stroemii + hediste_diversicolor + 
           phyllodocida + polychaeta + pectinaria_sp + nephtys_ciliata,
         saduria_entomon_tot = saduria_entomon
         )
#> mutate: new variable 'amphipoda_tot' (double) with 397 unique values and 0% NA
#>         new variable 'bivalvia_tot' (double) with 257 unique values and 0% NA
#>         new variable 'clupeidae_tot' (double) with 285 unique values and 0% NA
#>         new variable 'sprattus_sprattus_tot' (double) with 986 unique values and 0% NA
#>         new variable 'clupea_harengus_tot' (double) with 646 unique values and 0% NA
#>         new variable 'gadiformes_tot' (double) with 7 unique values and 0% NA
#>         new variable 'gobiidae_tot' (double) with 163 unique values and 0% NA
#>         new variable 'mysidae_tot' (double) with 358 unique values and 0% NA
#>         new variable 'non_bio_tot' (double) with 61 unique values and 0% NA
#>         new variable 'other_crustacea_tot' (double) with 285 unique values and 0% NA
#>         new variable 'other_tot' (double) with 299 unique values and 0% NA
#>         new variable 'other_pisces_tot' (double) with 253 unique values and 0% NA
#>         new variable 'platichthys_flesus_tot' (double) with 16 unique values and 0% NA
#>         new variable 'polychaeta_tot' (double) with 335 unique values and 0% NA
#>         new variable 'saduria_entomon_tot' (double) with 626 unique values and 0% NA

length(unique(d$Unique.pred.id))
#> [1] 9432
str(d_wide2) # Check why they do not match in length / nrows
#> tibble [8,311 × 128] (S3: tbl_df/tbl/data.frame)
#>  $ unique_pred_id            : chr [1:8311] "2007_4_101_2007_11_7_101_DAN2_40G5_1_M_490_1028_1_COD" "2007_4_101_2007_11_7_101_DAN2_40G5_118_F_210_82_1_COD" "2007_4_101_2007_11_7_101_DAN2_40G5_45_F_480_784_1_COD" "2007_4_101_2007_11_7_101_DAN2_40G5_65_F_920_7140_1_COD" ...
#>  $ sprattus_sprattus         : num [1:8311] 5.49 0 0 26.59 0 ...
#>  $ bylgides_sarsi            : num [1:8311] 0 0.24 0 0 0 0 0 0 0 0 ...
#>  $ saduria_entomon           : num [1:8311] 0 0 0.9 0 0 0 0 0 0 0 ...
#>  $ clupea_harengus           : num [1:8311] 0 0 0 179.7 88.5 ...
#>  $ pisces                    : num [1:8311] 0 0 0 20.8 0 ...
#>  $ gadus_morhua              : num [1:8311] 0 0 0 0 114 ...
#>  $ cumacea                   : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ limecola_balthica         : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ mysis_mixta               : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ corophium_volutator       : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ halicryptus_spinulosus    : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ stone                     : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ pontoporeia_femorata      : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ enchelyopus_cimbrius      : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ clupeidae                 : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ gobiidae                  : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ neomysis_integer          : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ monoporeia_affinis        : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ gammarus_sp               : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ mysidae                   : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ scoloplos_armiger         : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ plastic                   : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ priapulus_caudatus        : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ bivalvia                  : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ hyperia_galba             : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ mytilus_sp                : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ crangon_crangon           : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ idotea_balthica           : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ trachurus_trachurus       : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ gasterosteus_aculeatus    : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ annelida                  : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ algae                     : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ platichthys_flesus        : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ priapulidae               : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ waste                     : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ praunus_flexuosus         : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ unidentified_mass         : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ merlangius_merlangus      : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ scales                    : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ gadidae                   : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ crustacea                 : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ amphipoda                 : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ spine                     : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ pleuronectes_platessa     : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ anguilla_anguilla         : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ empty                     : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ mollusca                  : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ pungitius_pungitius       : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ diastylis_rathkei         : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ terebellides_stroemii     : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ palaemon_sp               : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ ammodytes_tobianus        : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ na                        : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ cottidae                  : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ hediste_diversicolor      : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ palaemon_elegans          : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ spinachia_spinachia       : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ zoarces_viviparus         : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ ammodytidae               : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ caridea                   : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ amphibalanus_improvisus   : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ myoxocephalus_quadricornis: num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ phyllodocida              : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ remains                   : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ palaemonidae              : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ carcinus_maenas           : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ gastropoda                : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ sand                      : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ cerastoderma_glaucum      : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ hyperoplus_lanceolatus    : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ mya_arenaria              : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ hydrobia_sp               : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ wood                      : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ pleuronectiformes         : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ scophthalmus_maximus      : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ polychaeta                : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ priapulida                : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ halicryptus               : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ neogobius_melanostomus    : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ gobius_niger              : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ digestive_tract           : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ pleuronectidae            : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ sprattus_sprattus_2       : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ gasterosteidae            : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ litter                    : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ mysida                    : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ carbon                    : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ copepoda                  : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ calanoida                 : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ belone_belone             : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ pontoporeiidae            : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ mucus                     : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ pectinaria_sp             : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ macoma_balthica           : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ remains_2                 : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ stone_2                   : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ carbon_2                  : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ nephtys_ciliata           : num [1:8311] 0 0 0 0 0 0 0 0 0 0 ...
#>   [list output truncated]

# Select only columns aggregated columns
d_wide3 <- d_wide2 %>% dplyr::select(c(1, 114:128))

colnames(d_wide3)  
#>  [1] "unique_pred_id"         "amphipoda_tot"          "bivalvia_tot"          
#>  [4] "clupeidae_tot"          "sprattus_sprattus_tot"  "clupea_harengus_tot"   
#>  [7] "gadiformes_tot"         "gobiidae_tot"           "mysidae_tot"           
#> [10] "non_bio_tot"            "other_crustacea_tot"    "other_tot"             
#> [13] "other_pisces_tot"       "platichthys_flesus_tot" "polychaeta_tot"        
#> [16] "saduria_entomon_tot"

# Add back in other information 
d_sub <- d %>%
  dplyr::select(Year, Quarter, Cruise, Predator.code, X, Y, Lat, Long, Ices.rect, Pred.size.mm,
                Pred.weight.g, Unique.pred.id, source) %>% 
  rename("Predator.spec" = "Predator.code") %>% 
  distinct(Unique.pred.id, .keep_all = TRUE) %>%
  janitor::clean_names()  
#> rename: renamed one variable (Predator.spec)
#> distinct: removed 12,299 rows (57%), 9,432 rows remaining

d_wide4 <- left_join(d_wide3, d_sub) # Why missing 1,121 IDs?
#> Joining, by = "unique_pred_id"
#> left_join: added 12 columns (year, quarter, cruise, predator_spec, x, …)
#>            > rows only in x       0
#>            > rows only in y  (1,121)
#>            > matched rows     8,311
#>            >                 =======
#>            > rows total       8,311

#d_wide4 %>% group_by(unique_pred_id) %>% mutate(n = n()) %>% ungroup() %>% distinct(n)
colnames(d_wide4)
#>  [1] "unique_pred_id"         "amphipoda_tot"          "bivalvia_tot"          
#>  [4] "clupeidae_tot"          "sprattus_sprattus_tot"  "clupea_harengus_tot"   
#>  [7] "gadiformes_tot"         "gobiidae_tot"           "mysidae_tot"           
#> [10] "non_bio_tot"            "other_crustacea_tot"    "other_tot"             
#> [13] "other_pisces_tot"       "platichthys_flesus_tot" "polychaeta_tot"        
#> [16] "saduria_entomon_tot"    "year"                   "quarter"               
#> [19] "cruise"                 "predator_spec"          "x"                     
#> [22] "y"                      "lat"                    "long"                  
#> [25] "ices_rect"              "pred_size_mm"           "pred_weight_g"         
#> [28] "source"
head(data.frame(d_wide4))
#>                                           unique_pred_id amphipoda_tot
#> 1  2007_4_101_2007_11_7_101_DAN2_40G5_1_M_490_1028_1_COD          0.00
#> 2  2007_4_101_2007_11_7_101_DAN2_40G5_118_F_210_82_1_COD          0.24
#> 3  2007_4_101_2007_11_7_101_DAN2_40G5_45_F_480_784_1_COD          0.00
#> 4 2007_4_101_2007_11_7_101_DAN2_40G5_65_F_920_7140_1_COD          0.00
#> 5 2007_4_101_2007_11_7_101_DAN2_40G5_66_M_740_3812_1_COD          0.00
#> 6 2007_4_101_2007_11_7_101_DAN2_40G5_70_F_590_2198_1_COD          0.00
#>   bivalvia_tot clupeidae_tot sprattus_sprattus_tot clupea_harengus_tot
#> 1            0             0                  5.49                0.00
#> 2            0             0                  0.00                0.00
#> 3            0             0                  0.00                0.00
#> 4            0             0                 26.59              179.66
#> 5            0             0                  0.00               88.51
#> 6            0             0                 27.32              163.47
#>   gadiformes_tot gobiidae_tot mysidae_tot non_bio_tot other_crustacea_tot
#> 1              0            0           0           0                   0
#> 2              0            0           0           0                   0
#> 3              0            0           0           0                   0
#> 4              0            0           0           0                   0
#> 5              0            0           0           0                   0
#> 6              0            0           0           0                   0
#>   other_tot other_pisces_tot platichthys_flesus_tot polychaeta_tot
#> 1         0                0                      0           0.00
#> 2         0                0                      0           0.24
#> 3         0                0                      0           0.00
#> 4         0                0                      0           0.00
#> 5         0                0                      0           0.00
#> 6         0                0                      0           0.00
#>   saduria_entomon_tot year quarter cruise predator_spec        x        y
#> 1                 0.0 2007       4 BITS-2           COD 520.8793 6176.115
#> 2                 0.0 2007       4 BITS-2           COD 520.8793 6176.115
#> 3                 0.9 2007       4 BITS-2           COD 520.8793 6176.115
#> 4                 0.0 2007       4 BITS-2           COD 520.8793 6176.115
#> 5                 0.0 2007       4 BITS-2           COD 520.8793 6176.115
#> 6                 0.0 2007       4 BITS-2           COD 520.8793 6176.115
#>        lat     long ices_rect pred_size_mm pred_weight_g      source
#> 1 55.73032 15.33247      40G5          490          1028 Kevin's MSC
#> 2 55.73032 15.33247      40G5          210            82 Kevin's MSC
#> 3 55.73032 15.33247      40G5          480           784 Kevin's MSC
#> 4 55.73032 15.33247      40G5          920          7140 Kevin's MSC
#> 5 55.73032 15.33247      40G5          740          3812 Kevin's MSC
#> 6 55.73032 15.33247      40G5          590          2198 Kevin's MSC

colnames(d_wide4)
#>  [1] "unique_pred_id"         "amphipoda_tot"          "bivalvia_tot"          
#>  [4] "clupeidae_tot"          "sprattus_sprattus_tot"  "clupea_harengus_tot"   
#>  [7] "gadiformes_tot"         "gobiidae_tot"           "mysidae_tot"           
#> [10] "non_bio_tot"            "other_crustacea_tot"    "other_tot"             
#> [13] "other_pisces_tot"       "platichthys_flesus_tot" "polychaeta_tot"        
#> [16] "saduria_entomon_tot"    "year"                   "quarter"               
#> [19] "cruise"                 "predator_spec"          "x"                     
#> [22] "y"                      "lat"                    "long"                  
#> [25] "ices_rect"              "pred_size_mm"           "pred_weight_g"         
#> [28] "source"
d_wide4[1, 2]
#> # A tibble: 1 x 1
#>   amphipoda_tot
#>           <dbl>
#> 1             0
d_wide4[1, 16]
#> # A tibble: 1 x 1
#>   saduria_entomon_tot
#>                 <dbl>
#> 1                   0

d_wide4 <- d_wide4 %>% mutate(tot_prey_biom = rowSums(.[2:16]))
#> mutate: new variable 'tot_prey_biom' (double) with 2,467 unique values and 0% NA

# Now, group by survey and Quarter, select only the main groups (as Haase) and see if it resemlbes the
# previous results.

fle <- d_wide4 %>%
  filter(predator_spec == "FLE") %>% 
  mutate(pred_cm = pred_size_mm/10,
         pred_cm_class = cut(pred_cm, breaks = c(0, 9, 20, 30, 200), right = TRUE))
#> filter: removed 6,054 rows (73%), 2,257 rows remaining
#> mutate: new variable 'pred_cm' (double) with 113 unique values and <1% NA
#>         new variable 'pred_cm_class' (factor) with 5 unique values and <1% NA

head(data.frame(fle), 20)
#>         unique_pred_id amphipoda_tot bivalvia_tot clupeidae_tot
#> 1    2015_2_1010_4_FLE             0         2.60             0
#> 2  2015_2_1010_402_FLE             0         0.00             0
#> 3     2015_2_108_3_FLE             0         0.00             0
#> 4   2015_2_202_403_FLE             0         5.46             0
#> 5     2015_2_203_1_FLE             0         1.74             0
#> 6     2015_2_204_2_FLE             0         0.04             0
#> 7    2015_2_501_10_FLE             0         0.99             0
#> 8    2015_2_501_60_FLE             0         0.00             0
#> 9     2015_2_501_8_FLE             0         0.00             0
#> 10    2015_2_501_9_FLE             0         0.00             0
#> 11    2015_2_504_5_FLE             0         0.00             0
#> 12   2015_2_504_58_FLE             0         0.00             0
#> 13    2015_2_504_6_FLE             0         0.00             0
#> 14    2015_2_504_7_FLE             0         0.00             0
#> 15   2015_2_506_11_FLE             0         0.00             0
#> 16   2015_2_506_12_FLE             0         0.00             0
#> 17   2015_2_506_13_FLE             0         0.00             0
#> 18   2015_2_506_14_FLE             0         0.00             0
#> 19   2015_2_506_15_FLE             0         0.00             0
#> 20   2015_2_506_16_FLE             0         0.00             0
#>    sprattus_sprattus_tot clupea_harengus_tot gadiformes_tot gobiidae_tot
#> 1                      0                   0              0            0
#> 2                      0                   0              0            0
#> 3                      0                   0              0            0
#> 4                      0                   0              0            0
#> 5                      0                   0              0            0
#> 6                      0                   0              0            0
#> 7                      0                   0              0            0
#> 8                      0                   0              0            0
#> 9                      0                   0              0            0
#> 10                     0                   0              0            0
#> 11                     0                   0              0            0
#> 12                     0                   0              0            0
#> 13                     0                   0              0            0
#> 14                     0                   0              0            0
#> 15                     0                   0              0            0
#> 16                     0                   0              0            0
#> 17                     0                   0              0            0
#> 18                     0                   0              0            0
#> 19                     0                   0              0            0
#> 20                     0                   0              0            0
#>    mysidae_tot non_bio_tot other_crustacea_tot other_tot other_pisces_tot
#> 1            0        0.19                0.00      0.00             2.93
#> 2            0        0.00                0.00      0.00             0.00
#> 3            0        0.00                0.00      0.00             0.00
#> 4            0        0.13                0.00      0.00             0.00
#> 5            0        0.00                0.00      0.00             0.00
#> 6            0        0.00                0.00      0.00             0.00
#> 7            0        0.00                0.00      0.00             0.00
#> 8            0        0.00                0.00      0.00             0.00
#> 9            0        0.00                0.01      0.03             0.00
#> 10           0        0.00                0.00      0.00             0.00
#> 11           0        0.00                0.02      0.00             0.00
#> 12           0        0.00                0.00      0.02             0.00
#> 13           0        0.00                0.00      0.00             0.00
#> 14           0        0.00                0.00      0.00             0.00
#> 15           0        0.00                0.00      0.00             0.00
#> 16           0        0.00                0.10      0.00             0.00
#> 17           0        0.00                0.00      0.00             0.00
#> 18           0        0.00                0.00      0.00             0.00
#> 19           0        0.00                0.00      0.00             0.00
#> 20           0        0.00                0.00      0.00             0.00
#>    platichthys_flesus_tot polychaeta_tot saduria_entomon_tot year quarter
#> 1                       0              0                0.00 2015       2
#> 2                       0              0                0.00 2015       2
#> 3                       0              0                5.61 2015       2
#> 4                       0              0                1.02 2015       2
#> 5                       0              0                3.10 2015       2
#> 6                       0              0                2.98 2015       2
#> 7                       0              0                0.00 2015       2
#> 8                       0              0                0.00 2015       2
#> 9                       0              0                0.00 2015       2
#> 10                      0              0                0.00 2015       2
#> 11                      0              0                0.08 2015       2
#> 12                      0              0                0.00 2015       2
#> 13                      0              0                0.00 2015       2
#> 14                      0              0                0.00 2015       2
#> 15                      0              0                0.00 2015       2
#> 16                      0              0                0.00 2015       2
#> 17                      0              0                0.00 2015       2
#> 18                      0              0                0.00 2015       2
#> 19                      0              0                0.00 2015       2
#> 20                      0              0                0.00 2015       2
#>     cruise predator_spec        x        y   lat  long ices_rect pred_size_mm
#> 1  INSPIRE           FLE 486.2897 6209.440 56.03 14.78      41G4          270
#> 2  INSPIRE           FLE 486.2897 6209.440 56.03 14.78      41G4          304
#> 3  INSPIRE           FLE 483.1607 6206.112 56.00 14.73      41G4          350
#> 4  INSPIRE           FLE 485.6369 6200.539 55.95 14.77      40G4          302
#> 5  INSPIRE           FLE 486.9129 6209.438 56.03 14.79      41G4          299
#> 6  INSPIRE           FLE 485.6517 6204.990 55.99 14.77      40G4          346
#> 7  INSPIRE           FLE 478.0250 6177.198 55.74 14.65      40G4          314
#> 8  INSPIRE           FLE 478.0250 6177.198 55.74 14.65      40G4          330
#> 9  INSPIRE           FLE 478.0250 6177.198 55.74 14.65      40G4          202
#> 10 INSPIRE           FLE 478.0250 6177.198 55.74 14.65      40G4          298
#> 11 INSPIRE           FLE 477.9969 6171.634 55.69 14.65      40G4          352
#> 12 INSPIRE           FLE 477.9969 6171.634 55.69 14.65      40G4          313
#> 13 INSPIRE           FLE 477.9969 6171.634 55.69 14.65      40G4          328
#> 14 INSPIRE           FLE 477.9969 6171.634 55.69 14.65      40G4          250
#> 15 INSPIRE           FLE 479.2754 6176.079 55.73 14.67      40G4          299
#> 16 INSPIRE           FLE 479.2754 6176.079 55.73 14.67      40G4          280
#> 17 INSPIRE           FLE 479.2754 6176.079 55.73 14.67      40G4          349
#> 18 INSPIRE           FLE 479.2754 6176.079 55.73 14.67      40G4          335
#> 19 INSPIRE           FLE 479.2754 6176.079 55.73 14.67      40G4          340
#> 20 INSPIRE           FLE 479.2754 6176.079 55.73 14.67      40G4          337
#>    pred_weight_g      source tot_prey_biom pred_cm pred_cm_class
#> 1            236 Kevin's MSC          5.72    27.0       (20,30]
#> 2            476 Kevin's MSC          0.00    30.4      (30,200]
#> 3            334 Kevin's MSC          5.61    35.0      (30,200]
#> 4            264 Kevin's MSC          6.61    30.2      (30,200]
#> 5            246 Kevin's MSC          4.84    29.9       (20,30]
#> 6            346 Kevin's MSC          3.02    34.6      (30,200]
#> 7            305 Kevin's MSC          0.99    31.4      (30,200]
#> 8            548 Kevin's MSC          0.00    33.0      (30,200]
#> 9             89 Kevin's MSC          0.04    20.2       (20,30]
#> 10           251 Kevin's MSC          0.00    29.8       (20,30]
#> 11           325 Kevin's MSC          0.10    35.2      (30,200]
#> 12           319 Kevin's MSC          0.02    31.3      (30,200]
#> 13           367 Kevin's MSC          0.00    32.8      (30,200]
#> 14           172 Kevin's MSC          0.00    25.0       (20,30]
#> 15           196 Kevin's MSC          0.00    29.9       (20,30]
#> 16           184 Kevin's MSC          0.10    28.0       (20,30]
#> 17           389 Kevin's MSC          0.00    34.9      (30,200]
#> 18           368 Kevin's MSC          0.00    33.5      (30,200]
#> 19           371 Kevin's MSC          0.00    34.0      (30,200]
#> 20           380 Kevin's MSC          0.00    33.7      (30,200]
unique(fle$pred_cm_class)
#> [1] (20,30]  (30,200] (9,20]   (0,9]    <NA>    
#> Levels: (0,9] (9,20] (20,30] (30,200]

cod <- d_wide4 %>%
  filter(predator_spec == "COD") %>% 
  mutate(pred_cm = pred_size_mm/10,
         pred_cm_class = cut(pred_cm, breaks = c(0, 6, 20, 30, 40, 50, 200), right = TRUE))
#> filter: removed 2,257 rows (27%), 6,054 rows remaining
#> mutate: new variable 'pred_cm' (double) with 222 unique values and <1% NA
#>         new variable 'pred_cm_class' (factor) with 7 unique values and <1% NA

head(data.frame(cod), 20)
#>                                             unique_pred_id amphipoda_tot
#> 1    2007_4_101_2007_11_7_101_DAN2_40G5_1_M_490_1028_1_COD          0.00
#> 2    2007_4_101_2007_11_7_101_DAN2_40G5_118_F_210_82_1_COD          0.24
#> 3    2007_4_101_2007_11_7_101_DAN2_40G5_45_F_480_784_1_COD          0.00
#> 4   2007_4_101_2007_11_7_101_DAN2_40G5_65_F_920_7140_1_COD          0.00
#> 5   2007_4_101_2007_11_7_101_DAN2_40G5_66_M_740_3812_1_COD          0.00
#> 6   2007_4_101_2007_11_7_101_DAN2_40G5_70_F_590_2198_1_COD          0.00
#> 7   2007_4_101_2007_11_7_101_DAN2_40G5_72_F_540_1312_1_COD          0.00
#> 8   2007_4_101_2007_11_7_101_DAN2_40G5_73_F_560_1618_1_COD          0.00
#> 9   2007_4_101_2007_11_7_101_DAN2_40G5_75_F_780_5026_1_COD          0.00
#> 10  2007_4_101_2007_11_7_101_DAN2_40G5_76_M_590_1868_1_COD          0.00
#> 11  2007_4_101_2007_11_7_101_DAN2_40G5_77_F_590_1320_1_COD          0.00
#> 12  2007_4_101_2007_11_7_101_DAN2_40G5_81_M_630_2002_1_COD          0.00
#> 13  2007_4_101_2007_11_7_101_DAN2_40G5_82_F_680_2934_1_COD          0.00
#> 14  2007_4_101_2007_11_7_101_DAN2_40G5_84_F_580_1382_1_COD          0.00
#> 15  2007_4_101_2007_11_7_101_DAN2_40G5_86_F_580_1360_1_COD          0.00
#> 16  2007_4_101_2007_11_7_101_DAN2_40G5_87_F_540_1380_1_COD          0.00
#> 17  2007_4_101_2007_11_7_101_DAN2_40G5_94_M_540_1536_1_COD          0.00
#> 18  2007_4_101_2007_11_7_101_DAN2_40G5_95_M_590_1714_1_COD          0.00
#> 19   2007_4_103_2007_11_7_103_DAN2_40G5_1_M_740_3350_1_COD          0.00
#> 20 2007_4_103_2007_11_7_103_DAN2_40G5_122_F_650_2588_1_COD          0.00
#>    bivalvia_tot clupeidae_tot sprattus_sprattus_tot clupea_harengus_tot
#> 1             0             0                  5.49                0.00
#> 2             0             0                  0.00                0.00
#> 3             0             0                  0.00                0.00
#> 4             0             0                 26.59              179.66
#> 5             0             0                  0.00               88.51
#> 6             0             0                 27.32              163.47
#> 7             0             0                  4.93                0.00
#> 8             0             0                 14.67                0.00
#> 9             0             0                 51.84              264.95
#> 10            0             0                  1.72                0.00
#> 11            0             0                 18.24                0.00
#> 12            0             0                 12.32               51.80
#> 13            0             0                 31.46              169.98
#> 14            0             0                  0.00                0.00
#> 15            0             0                 19.95                0.00
#> 16            0             0                  0.00                0.00
#> 17            0             0                 12.68               59.26
#> 18            0             0                 24.77               75.62
#> 19            0             0                  0.00                0.00
#> 20            0             0                  0.00                0.00
#>    gadiformes_tot gobiidae_tot mysidae_tot non_bio_tot other_crustacea_tot
#> 1               0            0           0           0                   0
#> 2               0            0           0           0                   0
#> 3               0            0           0           0                   0
#> 4               0            0           0           0                   0
#> 5               0            0           0           0                   0
#> 6               0            0           0           0                   0
#> 7               0            0           0           0                   0
#> 8               0            0           0           0                   0
#> 9               0            0           0           0                   0
#> 10              0            0           0           0                   0
#> 11              0            0           0           0                   0
#> 12              0            0           0           0                   0
#> 13              0            0           0           0                   0
#> 14              0            0           0           0                   0
#> 15              0            0           0           0                   0
#> 16              0            0           0           0                   0
#> 17              0            0           0           0                   0
#> 18              0            0           0           0                   0
#> 19              0            0           0           0                   0
#> 20              0            0           0           0                   0
#>    other_tot other_pisces_tot platichthys_flesus_tot polychaeta_tot
#> 1          0                0                      0           0.00
#> 2          0                0                      0           0.24
#> 3          0                0                      0           0.00
#> 4          0                0                      0           0.00
#> 5          0                0                      0           0.00
#> 6          0                0                      0           0.00
#> 7          0                0                      0           0.00
#> 8          0                0                      0           0.00
#> 9          0                0                      0           0.00
#> 10         0                0                      0           0.00
#> 11         0                0                      0           0.00
#> 12         0                0                      0           0.00
#> 13         0                0                      0           0.00
#> 14         0                0                      0           0.00
#> 15         0                0                      0           0.00
#> 16         0                0                      0           0.00
#> 17         0                0                      0           0.00
#> 18         0                0                      0           0.00
#> 19         0                0                      0           0.00
#> 20         0                0                      0           0.00
#>    saduria_entomon_tot year quarter cruise predator_spec        x        y
#> 1                0.000 2007       4 BITS-2           COD 520.8793 6176.115
#> 2                0.000 2007       4 BITS-2           COD 520.8793 6176.115
#> 3                0.900 2007       4 BITS-2           COD 520.8793 6176.115
#> 4                0.000 2007       4 BITS-2           COD 520.8793 6176.115
#> 5                0.000 2007       4 BITS-2           COD 520.8793 6176.115
#> 6                0.000 2007       4 BITS-2           COD 520.8793 6176.115
#> 7                0.000 2007       4 BITS-2           COD 520.8793 6176.115
#> 8                0.000 2007       4 BITS-2           COD 520.8793 6176.115
#> 9                0.000 2007       4 BITS-2           COD 520.8793 6176.115
#> 10               0.000 2007       4 BITS-2           COD 520.8793 6176.115
#> 11               0.950 2007       4 BITS-2           COD 520.8793 6176.115
#> 12               0.000 2007       4 BITS-2           COD 520.8793 6176.115
#> 13               0.000 2007       4 BITS-2           COD 520.8793 6176.115
#> 14               0.000 2007       4 BITS-2           COD 520.8793 6176.115
#> 15               0.000 2007       4 BITS-2           COD 520.8793 6176.115
#> 16              14.590 2007       4 BITS-2           COD 520.8793 6176.115
#> 17               8.691 2007       4 BITS-2           COD 520.8793 6176.115
#> 18               0.000 2007       4 BITS-2           COD 520.8793 6176.115
#> 19              58.260 2007       4 BITS-2           COD 524.3767 6182.403
#> 20              25.120 2007       4 BITS-2           COD 524.3767 6182.403
#>         lat     long ices_rect pred_size_mm pred_weight_g      source
#> 1  55.73032 15.33247      40G5          490          1028 Kevin's MSC
#> 2  55.73032 15.33247      40G5          210            82 Kevin's MSC
#> 3  55.73032 15.33247      40G5          480           784 Kevin's MSC
#> 4  55.73032 15.33247      40G5          920          7140 Kevin's MSC
#> 5  55.73032 15.33247      40G5          740          3812 Kevin's MSC
#> 6  55.73032 15.33247      40G5          590          2198 Kevin's MSC
#> 7  55.73032 15.33247      40G5          540          1312 Kevin's MSC
#> 8  55.73032 15.33247      40G5          560          1618 Kevin's MSC
#> 9  55.73032 15.33247      40G5          780          5026 Kevin's MSC
#> 10 55.73032 15.33247      40G5          590          1868 Kevin's MSC
#> 11 55.73032 15.33247      40G5          590          1320 Kevin's MSC
#> 12 55.73032 15.33247      40G5          630          2002 Kevin's MSC
#> 13 55.73032 15.33247      40G5          680          2934 Kevin's MSC
#> 14 55.73032 15.33247      40G5          580          1382 Kevin's MSC
#> 15 55.73032 15.33247      40G5          580          1360 Kevin's MSC
#> 16 55.73032 15.33247      40G5          540          1380 Kevin's MSC
#> 17 55.73032 15.33247      40G5          540          1536 Kevin's MSC
#> 18 55.73032 15.33247      40G5          590          1714 Kevin's MSC
#> 19 55.78665 15.38872      40G5          740          3350 Kevin's MSC
#> 20 55.78665 15.38872      40G5          650          2588 Kevin's MSC
#>    tot_prey_biom pred_cm pred_cm_class
#> 1          5.490      49       (40,50]
#> 2          0.480      21       (20,30]
#> 3          0.900      48       (40,50]
#> 4        206.250      92      (50,200]
#> 5         88.510      74      (50,200]
#> 6        190.790      59      (50,200]
#> 7          4.930      54      (50,200]
#> 8         14.670      56      (50,200]
#> 9        316.790      78      (50,200]
#> 10         1.720      59      (50,200]
#> 11        19.190      59      (50,200]
#> 12        64.120      63      (50,200]
#> 13       201.440      68      (50,200]
#> 14         0.000      58      (50,200]
#> 15        19.950      58      (50,200]
#> 16        14.590      54      (50,200]
#> 17        80.631      54      (50,200]
#> 18       100.390      59      (50,200]
#> 19        58.260      74      (50,200]
#> 20        25.120      65      (50,200]
unique(cod$pred_cm_class)
#> [1] (40,50]  (20,30]  (50,200] (6,20]   (30,40]  (0,6]    <NA>    
#> Levels: (0,6] (6,20] (20,30] (30,40] (40,50] (50,200]

Reproduce Haase figures with biomass proportions

# Try to reproduce Haase et al!
# Aggregate across predator individuals and convert to long data for plotting
cod2 <- cod %>%
  drop_na(pred_cm_class) %>% # This is because not all predators have length
  filter(cruise %in% c("BITS", "BITS-1", "BITS-2")) %>%
  filter(year %in% c(2015, 2016, 2017)) %>%
  filter(source == "Kevin's MSC") %>% 
  filter(!pred_cm_class == "(0,6]") %>% 
  group_by(quarter, pred_cm_class) %>% 
  summarise(sprattus_sprattus_tot_all = sum(sprattus_sprattus_tot),
            saduria_entomon_tot_all = sum(saduria_entomon_tot),
            mysidae_tot_all = sum(mysidae_tot),
            clupea_harengus_tot_all = sum(clupea_harengus_tot),
            sprattus_sprattus_tot_all = sum(sprattus_sprattus_tot),
            gobiidae_tot_all = sum(gobiidae_tot),
            gadiformes_tot_all = sum(gadiformes_tot),
            other_pisces_tot_all = sum(other_pisces_tot),
            polychaeta_tot_all = sum(polychaeta_tot)) %>% 
  ungroup() %>% 
  pivot_longer(3:10) %>% 
  group_by(quarter, pred_cm_class) %>% 
  mutate(percent = value/sum(value)) %>%
  ungroup()
#> drop_na: removed 24 rows (<1%), 6,030 rows remaining
#> filter: removed 749 rows (12%), 5,281 rows remaining
#> filter: removed 3,158 rows (60%), 2,123 rows remaining
#> filter: removed 283 rows (13%), 1,840 rows remaining
#> filter: removed 16 rows (1%), 1,824 rows remaining
#> group_by: 2 grouping variables (quarter, pred_cm_class)
#> summarise: now 10 rows and 10 columns, one group variable remaining (quarter)
#> ungroup: no grouping variables
#> pivot_longer: reorganized (sprattus_sprattus_tot_all, saduria_entomon_tot_all, mysidae_tot_all, clupea_harengus_tot_all, gobiidae_tot_all, …) into (name, value) [was 10x10, now 80x4]
#> group_by: 2 grouping variables (quarter, pred_cm_class)
#> mutate (grouped): new variable 'percent' (double) with 65 unique values and 0% NA
#> ungroup: no grouping variables
  
# So much herring for small cod?
# A 20 cm cod weighs approx 100 g. Can it have 20g of herring in the stomach?
#small_cod <-
d %>%
  filter(Pred.size.mm < 200 & Prey.sp.code %in% c("clupea harengus", "Clupea harengus")) %>% 
  dplyr::select(Prey.weight, Prey.sp.code, Unique.pred.id, Pred.size.mm)
#> filter: removed 21,728 rows (>99%), 3 rows remaining
#> # A tibble: 3 x 4
#>   Prey.weight Prey.sp.code   Unique.pred.id                         Pred.size.mm
#>         <dbl> <chr>          <chr>                                         <dbl>
#> 1       20.0  Clupea hareng… 2016_1_73_735_COD                               180
#> 2        4.48 Clupea hareng… 2007_4_126_2007_11_8_126_DAN2_39G4_1_…          160
#> 3        9.18 Clupea hareng… 2007_4_126_2007_11_8_126_DAN2_39G4_4_…          170

d_wide4 %>% filter(unique_pred_id == "2016_1_73_735_COD") %>% as.data.frame()
#> filter: removed 8,310 rows (>99%), one row remaining
#>      unique_pred_id amphipoda_tot bivalvia_tot clupeidae_tot
#> 1 2016_1_73_735_COD          0.05            0             0
#>   sprattus_sprattus_tot clupea_harengus_tot gadiformes_tot gobiidae_tot
#> 1                     0               19.99              0            0
#>   mysidae_tot non_bio_tot other_crustacea_tot other_tot other_pisces_tot
#> 1           0           0                   0         0                0
#>   platichthys_flesus_tot polychaeta_tot saduria_entomon_tot year quarter cruise
#> 1                      0           0.05                   0 2016       1   BITS
#>   predator_spec        x        y    lat   long ices_rect pred_size_mm
#> 1           COD 613.9639 6271.507 56.574 16.855      42G6          180
#>   pred_weight_g      source tot_prey_biom
#> 1            46 Kevin's MSC         20.09
d %>% filter(Unique.pred.id == "2016_1_73_735_COD") %>% as.data.frame()
#> filter: removed 21,729 rows (>99%), 2 rows remaining
#>   SubDiv Year Month Day N.with.food N.regurgitated N.skeletal N.empty HN Sample
#> 1     27 2016     2  28           1             NA         NA      NA 73    735
#> 2     27 2016     2  28          NA             NA          1      NA 73    735
#>   Length.code    Prey.sp.code Prey.size Prey.weight Prey.nr Stage.digestion
#> 1          Lt Clupea harengus       150       19.99       1               1
#> 2        <NA>  Bylgides sarsi        NA        0.05      NA               2
#>   Comment Parasites.in.stomach Quarter Coun. Ship Cruise Pred.size.mm
#> 1    <NA>                 <NA>       1   SWE   NA   BITS          180
#> 2    <NA>                 <NA>       1   SWE   NA   BITS          180
#>   Pred.weight.g Predator.gutted.weight Predator.code Sex Maturity Gall.content
#> 1            46                     42           COD   M        2            4
#> 2            46                     42           COD   M        2            4
#>   Q.year Age Date        Index    Lat   Long Depth.catch Ices.rect      source
#> 1 1_2016   2 <NA> OXBH.2016.73 56.574 16.855          63      42G6 Kevin's MSC
#> 2 1_2016   2 <NA> OXBH.2016.73 56.574 16.855          63      42G6 Kevin's MSC
#>   Number Sample.type transect Depthstep Gonad.weight.roundfish
#> 1     NA          NA       NA        NA                     NA
#> 2     NA          NA       NA        NA                     NA
#>   Liver.weight.roundfish Stomach.content Perc.stomac.content comments Processed
#> 1                     NA              NA                  NA       NA        NA
#> 2                     NA              NA                  NA       NA        NA
#>      Unique.pred.id        X        Y
#> 1 2016_1_73_735_COD 613.9639 6271.507
#> 2 2016_1_73_735_COD 613.9639 6271.507

# Check sample sizes per size class (numbers below bars in Kevins figures)
cod %>%
  drop_na(pred_cm_class) %>% # This is because not all predators have length
  filter(cruise %in% c("BITS", "BITS-1", "BITS-2")) %>%
  filter(year %in% c(2015, 2016, 2017)) %>%
  filter(source == "Kevin's MSC") %>% 
  filter(!pred_cm_class == "(0,6]") %>% 
  group_by(quarter, pred_cm_class) %>% 
  summarise(n = n())
#> drop_na: removed 24 rows (<1%), 6,030 rows remaining
#> filter: removed 749 rows (12%), 5,281 rows remaining
#> filter: removed 3,158 rows (60%), 2,123 rows remaining
#> filter: removed 283 rows (13%), 1,840 rows remaining
#> filter: removed 16 rows (1%), 1,824 rows remaining
#> group_by: 2 grouping variables (quarter, pred_cm_class)
#> summarise: now 10 rows and 3 columns, one group variable remaining (quarter)
#> # A tibble: 10 x 3
#> # Groups:   quarter [2]
#>    quarter pred_cm_class     n
#>      <dbl> <fct>         <int>
#>  1       1 (6,20]          137
#>  2       1 (20,30]         275
#>  3       1 (30,40]         303
#>  4       1 (40,50]         194
#>  5       1 (50,200]        120
#>  6       4 (6,20]          128
#>  7       4 (20,30]         266
#>  8       4 (30,40]         236
#>  9       4 (40,50]         132
#> 10       4 (50,200]         33
# I seem to have more data...  

# Check Kevin's raw data... 
kev <- read.csv("/Users/maxlindmark/Desktop/R_STUDIO_PROJECTS/bentfish/data/stomach_data/Master_Kevin_updated.csv", sep = ";")
head(kev)
#>   X coun. ship cruise sample.typ      sd ices rect transect       date  qyear
#> 1 1   SWE DANS   BITS      trawl SD26-28 43G8 4363          23.11.2015 4_2015
#> 2 2   SWE DANS   BITS      trawl SD26-28 43G8 4363          23.11.2015 4_2015
#> 3 3   SWE DANS   BITS      trawl SD26-28 43G8 4363          23.11.2015 4_2015
#> 4 4   SWE DANS   BITS      trawl SD26-28 43G8 4363          23.11.2015 4_2015
#> 5 5   SWE DANS   BITS      trawl SD26-28 43G8 4363          23.11.2015 4_2015
#> 6 6   SWE DANS   BITS      trawl SD26-28 43G8 4363          23.11.2015 4_2015
#>   year month day quarter station.haul        index    lat   long depth
#> 1 2015    11  23       4           19 OXBH.2015.19 57,005 18,806    72
#> 2 2015    11  23       4           19 OXBH.2015.19 57,005 18,806    72
#> 3 2015    11  23       4           19 OXBH.2015.19 57,005 18,806    72
#> 4 2015    11  23       4           19 OXBH.2015.19 57,005 18,806    72
#> 5 2015    11  23       4           19 OXBH.2015.19 57,005 18,806    72
#> 6 2015    11  23       4           19 OXBH.2015.19 57,005 18,806    72
#>   depthstep sample pred.code pred.size length lengthclass pred.weight
#> 1      >30m    406       COD       190     19      Lc6-20          70
#> 2      >30m    385       COD       320     32     Lc31-40         246
#> 3      >30m    389       COD       330     33     Lc31-40         330
#> 4      >30m    396       COD       300     30     Lc21-30         193
#> 5      >30m    381       COD       310     31     Lc31-40         244
#> 6      >30m    400       COD       250     25     Lc21-30         136
#>   pred.gutted.weight age sex maturity gonad.weight liver.weight gall.bladder
#> 1                 60   1   M        2           NA            5            1
#> 2                213   3   M        8           NA            9            1
#> 3                499   3   M        2           NA           27            4
#> 4                169   2   M        3           NA           10            4
#> 5                220   2   M        3           NA           11            4
#> 6                115   1   M        3           NA            8            2
#>   specimen.note n.full n.regurgitated n.skeletal n.empty stomachcontent
#> 1                    1             NA         NA      NA           0,87
#> 2                    1             NA         NA      NA           1,78
#> 3                   NA              1         NA      NA              0
#> 4                    1             NA         NA      NA           0,42
#> 5                    1             NA         NA      NA           0,06
#> 6                    1             NA         NA      NA            0,8
#>   perc.stomachcontent length.code    prey.sp.code prey.size prey.weight prey.nr
#> 1                1,24                 Mysis mixta        NA        0,87      19
#> 2                0,72          Ls Clupea harengus        85        1,78       1
#> 3                   0                                    NA           0        
#> 4                0,22                 Mysis mixta        NA        0,42      10
#> 5                0,02                 Mysis mixta        NA        0,06       2
#> 6                0,59                 Mysis mixta        NA         0,8      19
#>   stage.digestion      comment parasites.in.stomach processed
#> 1               1                                    original
#> 2               1 without head                       original
#> 3              NA fresh scales                           sure
#> 4               1                                    original
#> 5               1                                    original
#> 6               1                                    original

kev$prey.weight <- as.numeric(gsub(",", ".", gsub("\\.", "", kev$prey.weight)))

kev %>% 
  filter(lengthclass == "Lc6-20" & pred.code == "COD" & year %in% c(2015, 2016, 2017) & quarter == 1) %>% 
  drop_na(prey.sp.code) %>%
  drop_na(prey.weight) %>%
  group_by(prey.sp.code) %>% 
  summarise(tot_weight_by_prey = sum(prey.weight)) %>% 
  arrange(desc(tot_weight_by_prey)) %>% 
  as.data.frame()
#> filter: removed 15,869 rows (99%), 223 rows remaining
#> drop_na: no rows removed
#> drop_na: no rows removed
#> group_by: one grouping variable (prey.sp.code)
#> summarise: now 23 rows and 2 columns, ungrouped
#>              prey.sp.code tot_weight_by_prey
#> 1         Clupea harengus              19.99
#> 2             Mysis mixta               6.52
#> 3          Bylgides sarsi               4.52
#> 4         Saduria entomon               3.75
#> 5                Gobiidae               3.53
#> 6               Clupeidae               1.29
#> 7       Diastylis rathkei               1.11
#> 8         Crangon crangon               0.58
#> 9              Priapulida               0.40
#> 10           Gammarus sp.               0.28
#> 11                Mysidae               0.22
#> 12 Halicryptus spinulosus               0.17
#> 13                Cumacea               0.15
#> 14      Limecola balthica               0.11
#> 15               Bivalvia               0.10
#> 16   Pontoporeia femorata               0.07
#> 17                Remains               0.06
#> 18              Amphipoda               0.01
#> 19     Monoporeia affinis               0.01
#> 20            Mytilus sp.               0.01
#> 21       Neomysis integer               0.01
#> 22                                      0.00
#> 23               Copepoda               0.00

# OK so probably Kevin made some additional filter, because I also have a larger sample size

Plot potential response variables… Unlike Haase, I want to work with individual stomachs

# Plot and calculate proportion empty stomachs
head(data.frame(cod))
#>                                           unique_pred_id amphipoda_tot
#> 1  2007_4_101_2007_11_7_101_DAN2_40G5_1_M_490_1028_1_COD          0.00
#> 2  2007_4_101_2007_11_7_101_DAN2_40G5_118_F_210_82_1_COD          0.24
#> 3  2007_4_101_2007_11_7_101_DAN2_40G5_45_F_480_784_1_COD          0.00
#> 4 2007_4_101_2007_11_7_101_DAN2_40G5_65_F_920_7140_1_COD          0.00
#> 5 2007_4_101_2007_11_7_101_DAN2_40G5_66_M_740_3812_1_COD          0.00
#> 6 2007_4_101_2007_11_7_101_DAN2_40G5_70_F_590_2198_1_COD          0.00
#>   bivalvia_tot clupeidae_tot sprattus_sprattus_tot clupea_harengus_tot
#> 1            0             0                  5.49                0.00
#> 2            0             0                  0.00                0.00
#> 3            0             0                  0.00                0.00
#> 4            0             0                 26.59              179.66
#> 5            0             0                  0.00               88.51
#> 6            0             0                 27.32              163.47
#>   gadiformes_tot gobiidae_tot mysidae_tot non_bio_tot other_crustacea_tot
#> 1              0            0           0           0                   0
#> 2              0            0           0           0                   0
#> 3              0            0           0           0                   0
#> 4              0            0           0           0                   0
#> 5              0            0           0           0                   0
#> 6              0            0           0           0                   0
#>   other_tot other_pisces_tot platichthys_flesus_tot polychaeta_tot
#> 1         0                0                      0           0.00
#> 2         0                0                      0           0.24
#> 3         0                0                      0           0.00
#> 4         0                0                      0           0.00
#> 5         0                0                      0           0.00
#> 6         0                0                      0           0.00
#>   saduria_entomon_tot year quarter cruise predator_spec        x        y
#> 1                 0.0 2007       4 BITS-2           COD 520.8793 6176.115
#> 2                 0.0 2007       4 BITS-2           COD 520.8793 6176.115
#> 3                 0.9 2007       4 BITS-2           COD 520.8793 6176.115
#> 4                 0.0 2007       4 BITS-2           COD 520.8793 6176.115
#> 5                 0.0 2007       4 BITS-2           COD 520.8793 6176.115
#> 6                 0.0 2007       4 BITS-2           COD 520.8793 6176.115
#>        lat     long ices_rect pred_size_mm pred_weight_g      source
#> 1 55.73032 15.33247      40G5          490          1028 Kevin's MSC
#> 2 55.73032 15.33247      40G5          210            82 Kevin's MSC
#> 3 55.73032 15.33247      40G5          480           784 Kevin's MSC
#> 4 55.73032 15.33247      40G5          920          7140 Kevin's MSC
#> 5 55.73032 15.33247      40G5          740          3812 Kevin's MSC
#> 6 55.73032 15.33247      40G5          590          2198 Kevin's MSC
#>   tot_prey_biom pred_cm pred_cm_class
#> 1          5.49      49       (40,50]
#> 2          0.48      21       (20,30]
#> 3          0.90      48       (40,50]
#> 4        206.25      92      (50,200]
#> 5         88.51      74      (50,200]
#> 6        190.79      59      (50,200]
cod %>% 
  mutate(empty_stomach = ifelse(tot_prey_biom == 0, "Y", "N")) %>% 
  ggplot(., aes(empty_stomach, fill = empty_stomach)) + 
  geom_bar()
#> mutate: new variable 'empty_stomach' (character) with 2 unique values and 0% NA


# Plot "stomach" fullness, i.e., the weight of prey in stomach (relative to predator weight)
# Can also be called Feeding Ratio, FR
t <- cod %>% drop_na(pred_weight_g)
#> drop_na: removed 1,255 rows (21%), 4,799 rows remaining
t <- cod %>% drop_na(pred_cm)
#> drop_na: removed 24 rows (<1%), 6,030 rows remaining

# Next find the most abundance prey groups for cod and flounder (to see if the density has any effect
# on the biomass of common prey in their stomachs)
# To do that, I need long format again, group by prey item and summarise
colnames(cod)
#>  [1] "unique_pred_id"         "amphipoda_tot"          "bivalvia_tot"          
#>  [4] "clupeidae_tot"          "sprattus_sprattus_tot"  "clupea_harengus_tot"   
#>  [7] "gadiformes_tot"         "gobiidae_tot"           "mysidae_tot"           
#> [10] "non_bio_tot"            "other_crustacea_tot"    "other_tot"             
#> [13] "other_pisces_tot"       "platichthys_flesus_tot" "polychaeta_tot"        
#> [16] "saduria_entomon_tot"    "year"                   "quarter"               
#> [19] "cruise"                 "predator_spec"          "x"                     
#> [22] "y"                      "lat"                    "long"                  
#> [25] "ices_rect"              "pred_size_mm"           "pred_weight_g"         
#> [28] "source"                 "tot_prey_biom"          "pred_cm"               
#> [31] "pred_cm_class"
cod_important_prey <- cod %>%
  pivot_longer(2:16) %>% 
  group_by(name) %>% 
  summarise(tot_prey = sum(value)) %>% 
  mutate(percent = round(tot_prey / sum(tot_prey), digits = 2)*100) %>%  # calculate also percent of total
  arrange(desc(tot_prey)) %>% 
  mutate(spec = "cod")
#> pivot_longer: reorganized (amphipoda_tot, bivalvia_tot, clupeidae_tot, sprattus_sprattus_tot, clupea_harengus_tot, …) into (name, value) [was 6054x31, now 90810x18]
#> group_by: one grouping variable (name)
#> summarise: now 15 rows and 2 columns, ungrouped
#> mutate: new variable 'percent' (double) with 7 unique values and 0% NA
#> mutate: new variable 'spec' (character) with one unique value and 0% NA

# Same for flounder
fle_important_prey <- fle %>%
  pivot_longer(2:16) %>% 
  group_by(name) %>% 
  summarise(tot_prey = sum(value)) %>% 
  mutate(percent = round(tot_prey / sum(tot_prey), digits = 2)*100) %>%  # calculate also percent of total
  arrange(desc(tot_prey)) %>% 
  mutate(spec = "fle")
#> pivot_longer: reorganized (amphipoda_tot, bivalvia_tot, clupeidae_tot, sprattus_sprattus_tot, clupea_harengus_tot, …) into (name, value) [was 2257x31, now 33855x18]
#> group_by: one grouping variable (name)
#> summarise: now 15 rows and 2 columns, ungrouped
#> mutate: new variable 'percent' (double) with 9 unique values and 0% NA
#> mutate: new variable 'spec' (character) with one unique value and 0% NA

plotdat <- bind_rows(cod_important_prey, fle_important_prey)

ggplot(plotdat, 
       aes(name, percent, percent, fill = spec)) +#aes(fct_reorder(name, percent, .desc = TRUE), percent, fill = spec)) +
  geom_bar(stat = "identity") + 
  theme_classic(base_size = 16) + 
  #coord_flip() +
  theme(axis.text.x = element_text(angle = 90)) +
  scale_fill_brewer(palette = "Set2") + 
  NULL


# Ok, so the species that both prey feed on (even though very little(!)) are:
common_prey <- c("amphipoda_tot", "clupea_harengus_tot", "clupeidae_tot", "other_crustacea_tot",
                 "other_pisces_tot", "polychaeta_tot", "saduria_entomon_tot", "sprattus_sprattus_tot")


# Total feeding ratio, proportion of saduria and proportion of common prey!
# First for cod
cod_fr <- cod %>% 
  mutate(pred_weight_g = replace_na(pred_weight_g, -9)) %>% 
  mutate(pred_weight_g = ifelse(pred_weight_g == -9, 0.01*pred_cm^3, pred_weight_g)) %>% 
  mutate(FR_tot = (tot_prey_biom)/(pred_weight_g - tot_prey_biom),
         FR_saduria = (saduria_entomon_tot)/(pred_weight_g - tot_prey_biom)) %>%
  # mutate(FR_tot = tot_prey_biom/(pred_weight_g - tot_prey_biom),
  #        FR_saduria = saduria_entomon_tot/(pred_weight_g - tot_prey_biom)) %>%
  mutate(prop_saduria = saduria_entomon_tot/tot_prey_biom,
         prop_common = (amphipoda_tot + clupea_harengus_tot + clupeidae_tot + other_crustacea_tot +
                 other_pisces_tot + polychaeta_tot + saduria_entomon_tot + sprattus_sprattus_tot) / tot_prey_biom,
         common_tot = (amphipoda_tot + clupea_harengus_tot + clupeidae_tot + other_crustacea_tot +
                 other_pisces_tot + polychaeta_tot + saduria_entomon_tot + sprattus_sprattus_tot)) %>%
  filter(quarter %in% c(1, 4)) %>% 
  filter(lat < 58 & lat > 55) %>% 
  #filter(FR_tot > 0 & FR_tot < 1) #%>% # strange id: 2007_4_72_2007_11_6_72_DAN2_40G4_6_F_250_134_1_COD
  filter(FR_tot > 0)
#> mutate: changed 1,255 values (21%) of 'pred_weight_g' (1255 fewer NA)
#> mutate: changed 1,255 values (21%) of 'pred_weight_g' (24 new NA)
#> mutate: new variable 'FR_tot' (double) with 4,821 unique values and <1% NA
#>         new variable 'FR_saduria' (double) with 1,229 unique values and <1% NA
#> mutate: new variable 'prop_saduria' (double) with 836 unique values and 12% NA
#>         new variable 'prop_common' (double) with 1,242 unique values and 12% NA
#>         new variable 'common_tot' (double) with 2,086 unique values and 0% NA
#> filter: removed 91 rows (2%), 5,963 rows remaining
#> filter: removed 20 rows (<1%), 5,943 rows remaining
#> filter: removed 722 rows (12%), 5,221 rows remaining
  # I guess in theory the total stomach content could weigh more than the fish without food
  # dplyr::select(unique_pred_id, FR_tot) %>%
  # arrange(desc(FR_tot)) %>%
  # as.data.frame()

# Now flounder
# First for cod
fle_fr <- fle %>% 
  mutate(pred_weight_g = replace_na(pred_weight_g, -9)) %>% 
  mutate(pred_weight_g = ifelse(pred_weight_g == -9, 0.01*pred_cm^3, pred_weight_g)) %>% 
  mutate(FR_tot = (tot_prey_biom)/(pred_weight_g - tot_prey_biom),
         FR_saduria = (saduria_entomon_tot)/(pred_weight_g - tot_prey_biom)) %>%
  # mutate(FR_tot = tot_prey_biom/(pred_weight_g - tot_prey_biom),
  #        FR_saduria = saduria_entomon_tot/(pred_weight_g - tot_prey_biom)) %>%
  mutate(prop_saduria = saduria_entomon_tot/tot_prey_biom,
         prop_common = (amphipoda_tot + clupea_harengus_tot + clupeidae_tot + other_crustacea_tot +
                 other_pisces_tot + polychaeta_tot + saduria_entomon_tot + sprattus_sprattus_tot) / tot_prey_biom,
         common_tot = (amphipoda_tot + clupea_harengus_tot + clupeidae_tot + other_crustacea_tot +
                 other_pisces_tot + polychaeta_tot + saduria_entomon_tot + sprattus_sprattus_tot)) %>%
  filter(quarter %in% c(1, 4)) %>% 
  filter(lat < 58 & lat > 55) %>% 
  #filter(FR_tot > 0 & FR_tot < 1) #%>% # strange id: 2007_4_72_2007_11_6_72_DAN2_40G4_6_F_250_134_1_COD
  filter(FR_tot > 0)
#> mutate: changed 918 values (41%) of 'pred_weight_g' (918 fewer NA)
#> mutate: changed 918 values (41%) of 'pred_weight_g' (7 new NA)
#> mutate: new variable 'FR_tot' (double) with 1,601 unique values and <1% NA
#>         new variable 'FR_saduria' (double) with 640 unique values and <1% NA
#> mutate: new variable 'prop_saduria' (double) with 396 unique values and 25% NA
#>         new variable 'prop_common' (double) with 694 unique values and 25% NA
#>         new variable 'common_tot' (double) with 495 unique values and 0% NA
#> filter: removed 58 rows (3%), 2,199 rows remaining
#> filter: removed 10 rows (<1%), 2,189 rows remaining
#> filter: removed 522 rows (24%), 1,667 rows remaining
  # I guess in theory the total stomach content could weigh more than the fish without food
  # dplyr::select(unique_pred_id, FR_tot) %>%
  # arrange(desc(FR_tot)) %>%
  # as.data.frame()

# What is the mean FR per year
cod %>% 
  mutate(pred_weight_g = replace_na(pred_weight_g, -9)) %>% 
  mutate(pred_weight_g = ifelse(pred_weight_g == -9, 0.01*pred_cm^3, pred_weight_g)) %>% 
  mutate(FR_tot = (100 * tot_prey_biom)/(pred_weight_g - tot_prey_biom),
         FR_saduria = (100 * saduria_entomon_tot)/(pred_weight_g - tot_prey_biom)) %>%
  filter(FR_tot > -1) %>% 
  group_by(year) %>% 
  summarise(mean_FR_tot = mean(FR_tot)) %>% 
  ggplot(., aes(year, mean_FR_tot)) + 
  geom_point(size = 4) + 
  stat_smooth() + 
  theme_classic(base_size = 18)
#> mutate: changed 1,255 values (21%) of 'pred_weight_g' (1255 fewer NA)
#> mutate: changed 1,255 values (21%) of 'pred_weight_g' (24 new NA)
#> mutate: new variable 'FR_tot' (double) with 4,809 unique values and <1% NA
#>         new variable 'FR_saduria' (double) with 1,228 unique values and <1% NA
#> filter: removed 25 rows (<1%), 6,029 rows remaining
#> group_by: one grouping variable (year)
#> summarise: now 13 rows and 2 columns, ungrouped
#> `geom_smooth()` using method = 'loess' and formula 'y ~ x'


# Plot the distribution of FR
cod_fr %>% 
  ggplot(., aes(FR_tot)) + 
  geom_histogram() +  
  theme_classic(base_size = 14) +
  coord_cartesian(expand = 0) + 
  NULL
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.


plot_map_raster +
  geom_point(data = cod_fr, aes(x = x * 1000, y = y * 1000, color = log(FR_tot))) +
  scale_color_viridis() + 
  facet_wrap(cruise ~ quarter) +
  theme_plot()


# Instead plot the proportion of saduria
plot_map_raster +
  geom_point(data = cod_fr, aes(x = x * 1000, y = y * 1000, color = prop_saduria)) +
  scale_color_viridis() + 
  #facet_wrap(cruise ~ quarter) +
  theme_plot()


# Plot the distribution of proportions - color if the prop is zero!
cod_fr %>% 
  mutate(saduria_category = ifelse(prop_saduria == 0, "No", "Some")) %>%
  mutate(saduria_category = ifelse(prop_saduria == 1, "Only", saduria_category)) %>% 
  ggplot(., aes(prop_saduria, fill = saduria_category)) + 
  geom_histogram() +  
  scale_fill_brewer(palette = "Set2") + 
  theme_classic(base_size = 14) +
  coord_cartesian(expand = 0) + 
  annotate("text", label = paste("Mean proportion Saduria = ",
                                 round(mean(cod_fr$prop_saduria), digits = 3)), size = 6,
           x = 0.5, y = 2000) + # add the mean proportion as text!
  NULL
#> mutate: new variable 'saduria_category' (character) with 2 unique values and 0% NA
#> mutate: changed 273 values (5%) of 'saduria_category' (0 new NA)
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.


cod_fr %>% 
  mutate(saduria_category = ifelse(prop_saduria == 0, "No", "Some")) %>%
  mutate(saduria_category = ifelse(prop_saduria == 1, "Only", saduria_category)) %>% 
  ggplot(., aes(prop_saduria, fill = saduria_category)) + 
  geom_histogram() +  
  facet_wrap(~ year, scales = "free") + 
  scale_fill_brewer(palette = "Set2") + 
  theme_classic(base_size = 14) +
  coord_cartesian(expand = 0) + 
  NULL
#> mutate: new variable 'saduria_category' (character) with 2 unique values and 0% NA
#> mutate: changed 273 values (5%) of 'saduria_category' (0 new NA)
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Conduct test analyses

How does the proportion and total weight of prey per predator weight (prey here being total prey biomass, saduria biomass and biomass of common prey) relate to flounder, cod and total density of the two species? Do the same for flounder!

# First read the density models and predict the values
qwraps2::lazyload_cache_dir(path = "R/analysis/density_spatial_trend_models_cache/html")
#> Lazyloading: R/analysis/density_spatial_trend_models_cache/html/condition model with spatial trend
#> Lazyloading: R/analysis/density_spatial_trend_models_cache/html/fit cod covariates and field
#> Lazyloading: R/analysis/density_spatial_trend_models_cache/html/fit cod depth covariate
#> Lazyloading: R/analysis/density_spatial_trend_models_cache/html/fit cod no covars
#> Lazyloading: R/analysis/density_spatial_trend_models_cache/html/fit fle covariates and field
#> Lazyloading: R/analysis/density_spatial_trend_models_cache/html/fit fle depth covariate
#> Lazyloading: R/analysis/density_spatial_trend_models_cache/html/fit fle no covars
#> Lazyloading: R/analysis/density_spatial_trend_models_cache/html/make barrier spde mesh
#> Lazyloading: R/analysis/density_spatial_trend_models_cache/html/predict on grid
#> Lazyloading: R/analysis/density_spatial_trend_models_cache/html/spatial trends

# Make predictions onto the stomach data defined earlier (1 row per stomach and prey?)

# Add depth and rename coordinates to match the data used for fitting the density model
# Read the tifs
west <- raster("/Users/maxlindmark/Desktop/R_STUDIO_PROJECTS/cod_condition/data/depth_geo_tif/D5_2018_rgb-1.tif")
#plot(west)

east <- raster("/Users/maxlindmark/Desktop/R_STUDIO_PROJECTS/cod_condition/data/depth_geo_tif/D6_2018_rgb-1.tif")
# plot(east)

dep_rast <- raster::merge(west, east)

cod_fr$depth <- extract(dep_rast, cod_fr[, 24:23])
fle_fr$depth <- extract(dep_rast, fle_fr[, 24:23])

# Convert to depth (instead of elevation)
ggplot(cod_fr, aes(depth)) + geom_histogram()
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

cod_fr$depth <- (cod_fr$depth - max(drop_na(cod_fr)$depth)) *-1
#> drop_na: no rows removed
ggplot(cod_fr, aes(depth)) + geom_histogram()
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.


fle_fr$depth <- (fle_fr$depth - max(drop_na(fle_fr)$depth)) *-1
#> drop_na: no rows removed

ggplot(cod_fr, aes(x, y, color = depth)) + 
  geom_point() + 
  scale_color_viridis()


cod_fr <- cod_fr %>% rename("X" = "x", "Y" = "y")
#> rename: renamed 2 variables (X, Y)
cod_fr <- cod_fr %>% mutate(year = as.integer(year))
#> mutate: converted 'year' from double to integer (0 new NA)
cod_fr <- cod_fr %>% mutate(year = as.integer(year))
#> mutate: no changes

fle_fr <- fle_fr %>% rename("X" = "x", "Y" = "y")
#> rename: renamed 2 variables (X, Y)
fle_fr <- fle_fr %>% mutate(year = as.integer(year))
#> mutate: converted 'year' from double to integer (0 new NA)
fle_fr <- fle_fr %>% mutate(year = as.integer(year))
#> mutate: no changes

# Standardize depth to the DENSITY data, not this data set
density_dat <- readr::read_csv("https://raw.githubusercontent.com/maxlindmark/cod_condition/master/data/for_analysis/mdat_cpue.csv")
#> 
#> ── Column specification ────────────────────────────────────────────────────────
#> cols(
#>   density = col_double(),
#>   year = col_double(),
#>   lat = col_double(),
#>   lon = col_double(),
#>   quarter = col_double(),
#>   haul.id = col_character(),
#>   IDx = col_character(),
#>   ices_rect = col_character(),
#>   SD = col_double(),
#>   density_fle = col_double(),
#>   depth = col_double(),
#>   oxy = col_double(),
#>   temp = col_double(),
#>   X = col_double(),
#>   Y = col_double()
#> )

cod_fr <- cod_fr %>% mutate(depth_sc = (depth - mean(density_dat$depth)) / sd(density_dat$depth))
#> mutate: new variable 'depth_sc' (double) with 98 unique values and 0% NA
fle_fr <- fle_fr %>% mutate(depth_sc = (depth - mean(density_dat$depth)) / sd(density_dat$depth))
#> mutate: new variable 'depth_sc' (double) with 77 unique values and 0% NA

cod_fr2 <- cod_fr %>% filter(year < 2020) # Remove this after I've refitted the models with 2020!
#> filter: removed 258 rows (5%), 4,963 rows remaining
fle_fr2 <- fle_fr %>% filter(year < 2020) # Remove this after I've refitted the models with 2020!
#> filter: removed 433 rows (26%), 1,234 rows remaining

# Now predict using the density models
cod_stomach_predcod_density <- predict(mcod2, newdata = dplyr::select(cod_fr2, depth_sc, year, X, Y))
#> Warning: Detected a `s()` smoother. Smoothers are penalized in sdmTMB as
#> of version 0.0.19, but used to be unpenalized.
#> You no longer need to specify `k`, specifing `k` now indicates
#> an upper limit on the basis functions, and the degree of wiggliness
#> is determined by the data.
cod_stomach_predfle_density <- predict(mfle2, newdata = dplyr::select(cod_fr2, depth_sc, year, X, Y))
#> Warning: Detected a `s()` smoother. Smoothers are penalized in sdmTMB as
#> of version 0.0.19, but used to be unpenalized.
#> You no longer need to specify `k`, specifing `k` now indicates
#> an upper limit on the basis functions, and the degree of wiggliness
#> is determined by the data.

fle_stomach_predcod_density <- predict(mcod2, newdata = dplyr::select(fle_fr2, depth_sc, year, X, Y))
#> Warning: Detected a `s()` smoother. Smoothers are penalized in sdmTMB as
#> of version 0.0.19, but used to be unpenalized.
#> You no longer need to specify `k`, specifing `k` now indicates
#> an upper limit on the basis functions, and the degree of wiggliness
#> is determined by the data.
fle_stomach_predfle_density <- predict(mfle2, newdata = dplyr::select(fle_fr2, depth_sc, year, X, Y))
#> Warning: Detected a `s()` smoother. Smoothers are penalized in sdmTMB as
#> of version 0.0.19, but used to be unpenalized.
#> You no longer need to specify `k`, specifing `k` now indicates
#> an upper limit on the basis functions, and the degree of wiggliness
#> is determined by the data.

# head(dplyr::select(cod_fr2, X, Y))
# head(predcod_density)
# head(predfle_density)

cod_fr2$predcod_density <- exp(cod_stomach_predcod_density$est)
cod_fr2$predfle_density <- exp(cod_stomach_predfle_density$est)

fle_fr2$predcod_density <- exp(fle_stomach_predcod_density$est)
fle_fr2$predfle_density <- exp(fle_stomach_predfle_density$est)

# Add scales predicted densities to the stomach data
cod_fr2 <- cod_fr2 %>%
  mutate(predcod_density_sc = (predcod_density - mean(predcod_density))/sd(predcod_density),
         predfle_density_sc = (predfle_density - mean(predfle_density))/sd(predfle_density))
#> mutate: new variable 'predcod_density_sc' (double) with 396 unique values and 0% NA
#>         new variable 'predfle_density_sc' (double) with 397 unique values and 0% NA

fle_fr2 <- fle_fr2 %>%
  mutate(predcod_density_sc = (predcod_density - mean(predcod_density))/sd(predcod_density),
         predfle_density_sc = (predfle_density - mean(predfle_density))/sd(predfle_density))
#> mutate: new variable 'predcod_density_sc' (double) with 168 unique values and 0% NA
#>         new variable 'predfle_density_sc' (double) with 168 unique values and 0% NA

# Now plot stomach summaries against predicted densities
head(cod_fr2)
#> # A tibble: 6 x 42
#>   unique_pred_id amphipoda_tot bivalvia_tot clupeidae_tot sprattus_spratt…
#>   <chr>                  <dbl>        <dbl>         <dbl>            <dbl>
#> 1 2007_4_101_20…          0               0             0             5.49
#> 2 2007_4_101_20…          0.24            0             0             0   
#> 3 2007_4_101_20…          0               0             0             0   
#> 4 2007_4_101_20…          0               0             0            26.6 
#> 5 2007_4_101_20…          0               0             0             0   
#> 6 2007_4_101_20…          0               0             0            27.3 
#> # … with 37 more variables: clupea_harengus_tot <dbl>, gadiformes_tot <dbl>,
#> #   gobiidae_tot <dbl>, mysidae_tot <dbl>, non_bio_tot <dbl>,
#> #   other_crustacea_tot <dbl>, other_tot <dbl>, other_pisces_tot <dbl>,
#> #   platichthys_flesus_tot <dbl>, polychaeta_tot <dbl>,
#> #   saduria_entomon_tot <dbl>, year <int>, quarter <dbl>, cruise <chr>,
#> #   predator_spec <chr>, X <dbl>, Y <dbl>, lat <dbl>, long <dbl>,
#> #   ices_rect <chr>, pred_size_mm <dbl>, pred_weight_g <dbl>, source <chr>,
#> #   tot_prey_biom <dbl>, pred_cm <dbl>, pred_cm_class <fct>, FR_tot <dbl>,
#> #   FR_saduria <dbl>, prop_saduria <dbl>, prop_common <dbl>, common_tot <dbl>,
#> #   depth <dbl>, depth_sc <dbl>, predcod_density <dbl>, predfle_density <dbl>,
#> #   predcod_density_sc <dbl>, predfle_density_sc <dbl>

Plot cod stomach against density

# Common prey
p1 <- ggplot(cod_fr2, aes(predcod_density_sc, prop_common)) + 
  geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) + 
  stat_smooth(size = 2) + 
  theme_classic(base_size = 13) +
  ggtitle("prop common prey vs cod density")

p2 <- ggplot(cod_fr2, aes(predcod_density_sc, common_tot/pred_weight_g)) + 
  geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) + 
  stat_smooth(size = 2) + 
  theme_classic(base_size = 13) +
  ggtitle("g/g common prey vs cod density")

p3 <- ggplot(cod_fr2, aes(predfle_density_sc, prop_common)) + 
  geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) + 
  stat_smooth(size = 2) + 
  theme_classic(base_size = 13) +
  ggtitle("prop common prey vs fle density")

p4 <- ggplot(cod_fr2, aes(predfle_density_sc, common_tot/pred_weight_g)) + 
  geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) + 
  stat_smooth(size = 2) + 
  theme_classic(base_size = 13) +
  ggtitle("g/g common prey vs fle density")

p5 <- ggplot(cod_fr2, aes(predfle_density_sc + predcod_density_sc, prop_common)) + 
  geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) + 
  stat_smooth(size = 2) + 
  theme_classic(base_size = 13) +
  ggtitle("prop common prey vs fle + cod density")

p6 <- ggplot(cod_fr2, aes(predfle_density_sc + predcod_density_sc, common_tot/pred_weight_g)) + 
  geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) + 
  stat_smooth(size = 2) + 
  theme_classic(base_size = 13) +
  ggtitle("g/g common prey vs fle + coddensity")

(p1|p2) / (p3|p4) / (p5|p6) + plot_annotation(title = "Cod stomachs", 
                                                  theme = theme(plot.title = element_text(size = 16)))
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'


# Saduria
p7 <- ggplot(cod_fr2, aes(predcod_density_sc, prop_saduria)) + 
  geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) + 
  stat_smooth(size = 2) + 
  theme_classic(base_size = 13) +
  ggtitle("prop saduria vs cod density")

p8 <- ggplot(cod_fr2, aes(predcod_density_sc, saduria_entomon_tot/pred_weight_g)) + 
  geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) + 
  stat_smooth(size = 2) + 
  theme_classic(base_size = 13) +
  ggtitle("g/g saduria vs cod density")

p9 <- ggplot(cod_fr2, aes(predfle_density_sc, prop_saduria)) + 
  geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) + 
  stat_smooth(size = 2) + 
  theme_classic(base_size = 13) +
  ggtitle("prop saduria vs fle density")

p10 <- ggplot(cod_fr2, aes(predfle_density_sc, saduria_entomon_tot/pred_weight_g)) + 
  geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) + 
  stat_smooth(size = 2) + 
  theme_classic(base_size = 13) +
  ggtitle("g/g saduria vs fle density")

p11 <- ggplot(cod_fr2, aes(predfle_density_sc + predcod_density_sc, prop_saduria)) + 
  geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) + 
  stat_smooth(size = 2) + 
  theme_classic(base_size = 13) +
  ggtitle("prop saduria vs fle + cod density")

p12 <- ggplot(cod_fr2, aes(predfle_density_sc + predcod_density_sc, saduria_entomon_tot/pred_weight_g)) + 
  geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) + 
  stat_smooth(size = 2) + 
  theme_classic(base_size = 13) +
  ggtitle("g/g saduria vs fle  + coddensity")

(p7|p8) / (p9|p10) / (p11|p12) + plot_annotation(title = "Cod stomachs", 
                                                  theme = theme(plot.title = element_text(size = 16)))
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'


# Total stomach content (FR)
p13 <- ggplot(cod_fr2, aes(predcod_density_sc, FR_tot)) + 
  geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) + 
  stat_smooth(size = 2) + 
  theme_classic(base_size = 13) +
  ggtitle("FR vs cod density")

p14 <- ggplot(cod_fr2, aes(predfle_density_sc, FR_tot)) + 
  geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) + 
  stat_smooth(size = 2) + 
  theme_classic(base_size = 13) +
  ggtitle("FR vs fle density")

p15 <- ggplot(cod_fr2, aes(predfle_density_sc + predcod_density, FR_tot)) + 
  geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) + 
  stat_smooth(size = 2) + 
  theme_classic(base_size = 13) +
  ggtitle("FR vs fle  + coddensity")

(p13/p14/p15) + plot_annotation(title = "Cod stomachs", theme = theme(plot.title = element_text(size = 16)))
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'


# Closer inspection of some of the plots...
ggplot(cod_fr2, aes(predfle_density_sc, prop_saduria)) + 
  geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) + 
  stat_smooth(size = 2, method = "lm") + 
  theme_classic(base_size = 13) +
  ggtitle("prop saduria vs fle density")
#> `geom_smooth()` using formula 'y ~ x'


summary(lm(cod_fr2$prop_saduria ~ cod_fr2$predfle_density))
#> 
#> Call:
#> lm(formula = cod_fr2$prop_saduria ~ cod_fr2$predfle_density)
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -0.1330 -0.1321 -0.1271 -0.1001  0.9968 
#> 
#> Coefficients:
#>                           Estimate Std. Error t value Pr(>|t|)    
#> (Intercept)              1.330e-01  4.423e-03  30.082  < 2e-16 ***
#> cod_fr2$predfle_density -1.882e-04  3.076e-05  -6.119 1.02e-09 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 0.2942 on 4961 degrees of freedom
#> Multiple R-squared:  0.00749,    Adjusted R-squared:  0.00729 
#> F-statistic: 37.44 on 1 and 4961 DF,  p-value: 1.015e-09

ggplot(cod_fr2, aes(predfle_density, tot_prey_biom/pred_weight_g)) + 
  geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) + 
  stat_smooth(size = 2) + 
  theme_classic(base_size = 13) +
  ggtitle("g/g tot. stomach vs fle density")
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'


ggplot(cod_fr2, aes(predfle_density + predcod_density, tot_prey_biom/pred_weight_g)) + 
  geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) + 
  stat_smooth(size = 2) + 
  theme_classic(base_size = 13) +
  ggtitle("g/g tot. stomach vs fle  + coddensity")
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'

Plot flounder stomach against density

# Common prey
p1 <- ggplot(fle_fr2, aes(predcod_density_sc, prop_common)) + 
  geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) + 
  stat_smooth(size = 2) + 
  theme_classic(base_size = 13) +
  ggtitle("prop common prey vs cod density")

p2 <- ggplot(fle_fr2, aes(predcod_density_sc, common_tot/pred_weight_g)) + 
  geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) + 
  stat_smooth(size = 2) + 
  theme_classic(base_size = 13) +
  ggtitle("g/g common prey vs cod density")

p3 <- ggplot(fle_fr2, aes(predfle_density_sc, prop_common)) + 
  geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) + 
  stat_smooth(size = 2) + 
  theme_classic(base_size = 13) +
  ggtitle("prop common prey vs fle density")

p4 <- ggplot(fle_fr2, aes(predfle_density_sc, common_tot/pred_weight_g)) + 
  geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) + 
  stat_smooth(size = 2) + 
  theme_classic(base_size = 13) +
  ggtitle("g/g common prey vs fle density")

p5 <- ggplot(fle_fr2, aes(predfle_density_sc + predcod_density_sc, prop_common)) + 
  geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) + 
  stat_smooth(size = 2) + 
  theme_classic(base_size = 13) +
  ggtitle("prop common prey vs fle + cod density")

p6 <- ggplot(fle_fr2, aes(predfle_density_sc + predcod_density_sc, common_tot/pred_weight_g)) + 
  geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) + 
  stat_smooth(size = 2) + 
  theme_classic(base_size = 13) +
  ggtitle("g/g common prey vs fle + coddensity")

(p1|p2) / (p3|p4) / (p5|p6) + plot_annotation(title = "Flounder stomachs", 
                                                  theme = theme(plot.title = element_text(size = 16)))
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'


# Saduria
p7 <- ggplot(fle_fr2, aes(predcod_density_sc, prop_saduria)) + 
  geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) + 
  stat_smooth(size = 2) + 
  theme_classic(base_size = 13) +
  ggtitle("prop saduria vs cod density")

p8 <- ggplot(fle_fr2, aes(predcod_density_sc, saduria_entomon_tot/pred_weight_g)) + 
  geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) + 
  stat_smooth(size = 2) + 
  theme_classic(base_size = 13) +
  ggtitle("g/g saduria vs cod density")

p9 <- ggplot(fle_fr2, aes(predfle_density_sc, prop_saduria)) + 
  geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) + 
  stat_smooth(size = 2) + 
  theme_classic(base_size = 13) +
  ggtitle("prop saduria vs fle density")

p10 <- ggplot(fle_fr2, aes(predfle_density_sc, saduria_entomon_tot/pred_weight_g)) + 
  geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) + 
  stat_smooth(size = 2) + 
  theme_classic(base_size = 13) +
  ggtitle("g/g saduria vs fle density")

p11 <- ggplot(fle_fr2, aes(predfle_density_sc + predcod_density_sc, prop_saduria)) + 
  geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) + 
  stat_smooth(size = 2) + 
  theme_classic(base_size = 13) +
  ggtitle("prop saduria vs fle + cod density")

p12 <- ggplot(fle_fr2, aes(predfle_density_sc + predcod_density_sc, saduria_entomon_tot/pred_weight_g)) + 
  geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) + 
  stat_smooth(size = 2) + 
  theme_classic(base_size = 13) +
  ggtitle("g/g saduria vs fle  + coddensity")

(p7|p8) / (p9|p10) / (p11|p12) + plot_annotation(title = "Flounder stomachs", 
                                                  theme = theme(plot.title = element_text(size = 16)))
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'


# FR
p13 <- ggplot(fle_fr2, aes(predcod_density_sc, FR_tot)) + 
  geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) + 
  stat_smooth(size = 2) + 
  theme_classic(base_size = 13) +
  ggtitle("FR vs cod density")

p14 <- ggplot(fle_fr2, aes(predfle_density_sc, FR_tot)) + 
  geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) + 
  stat_smooth(size = 2) + 
  theme_classic(base_size = 13) +
  ggtitle("FR vs fle density")

p15 <- ggplot(fle_fr2, aes(predfle_density_sc + predcod_density, FR_tot)) + 
  geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) + 
  stat_smooth(size = 2) + 
  theme_classic(base_size = 13) +
  ggtitle("FR vs fle + coddensity")

(p13/p14/p15) + plot_annotation(title = "Flounder stomachs", theme = theme(plot.title = element_text(size = 16)))
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'

More in depth plots…

# Does it seem like proportions of the diet changes with density, but not the total amount of food?
# Check saduria and flounder for instance.
pp7 <- ggplot(fle_fr2, aes(predcod_density_sc, prop_saduria)) + 
  geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) + 
  stat_smooth(size = 2, method = "lm") + 
  theme_classic(base_size = 13) +
  ggtitle("prop saduria vs cod density")

pp8 <- ggplot(filter(fle_fr2, (saduria_entomon_tot/pred_weight_g) < 0.1),
             aes(predcod_density_sc, saduria_entomon_tot/pred_weight_g)) + 
  geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) + 
  stat_smooth(size = 2, method = "lm") + 
  theme_classic(base_size = 13) +
  ggtitle("g/g saduria vs cod density")
#> filter: removed 2 rows (<1%), 1,232 rows remaining

pp13 <- ggplot(filter(fle_fr2, (saduria_entomon_tot/pred_weight_g) < 0.1),
              aes(predcod_density_sc, FR_tot)) + 
  geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) + 
  stat_smooth(size = 2, method = "lm") + 
  theme_classic(base_size = 13) +
  ggtitle("g/g tot. stomach vs cod density")
#> filter: removed 2 rows (<1%), 1,232 rows remaining

ppfle <- pp7/pp8/pp13 +
  plot_layout(ncol = 3) + 
  plot_annotation(title = "Flounder stomachs", theme = theme(plot.title = element_text(size = 16)))

# And that cod density is unaffected (or less at least)
p7 <- ggplot(filter(cod_fr2, (saduria_entomon_tot/pred_weight_g) < 0.02 & predfle_density_sc < 2),
             aes(predfle_density_sc, prop_saduria)) + # Crop it up a little!
  geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) + 
  stat_smooth(size = 2, method = "lm") + 
  theme_classic(base_size = 13) +
  ggtitle("prop saduria vs fle density")
#> filter: removed 171 rows (3%), 4,792 rows remaining

p8 <- ggplot(filter(cod_fr2, (saduria_entomon_tot/pred_weight_g) < 0.02 & predfle_density_sc < 2),
             aes(predfle_density_sc, saduria_entomon_tot/pred_weight_g)) + 
  geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) + 
  stat_smooth(size = 2, method = "lm") + 
  theme_classic(base_size = 13) +
  ggtitle("g/g saduria vs fle density")
#> filter: removed 171 rows (3%), 4,792 rows remaining

p13 <- ggplot(filter(cod_fr2, (saduria_entomon_tot/pred_weight_g) < 0.02 & predfle_density_sc < 2 & FR_tot < 0.5),
              aes(predfle_density_sc, FR_tot)) + 
  geom_point(size = 2, shape = 21, fill = "black", color = "white", alpha = 0.8) + 
  stat_smooth(size = 2, method = "lm") + 
  theme_classic(base_size = 13) +
  ggtitle("g/g tot. stomach vs fle density")
#> filter: removed 179 rows (4%), 4,784 rows remaining

ppcod <- p7/p8/p13 +
  plot_layout(ncol = 3) +
  plot_annotation(title = "Cod stomachs", theme = theme(plot.title = element_text(size = 16)))

ppfle / ppcod
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'


# What does it mean? Yes, the proportion *may* be negatively affected by densities, 
# but the food content is not, which I guess would matter for condition

Try and fit sdmTMB models with densities as covariates…

# Cod 
pcod_spde <- make_mesh(cod_fr2, c("X", "Y"), n_knots = 150, type = "kmeans", seed = 42)
plot(pcod_spde)


pfle_spde <- make_mesh(fle_fr2, c("X", "Y"), n_knots = 70, type = "kmeans", seed = 42)
plot(pfle_spde)


# (add small values if 0, subtract small values if 1) for Beta regressions
cod_fr2 <- cod_fr2 %>%
  mutate(saduria_entomon_per_mass = saduria_entomon_tot/pred_weight_g,
         prop_saduria2 = ifelse(prop_saduria == 0, prop_saduria + 0.0001, prop_saduria),
         prop_saduria2 = ifelse(prop_saduria == 1, prop_saduria - 0.0001, prop_saduria2))
#> mutate: new variable 'saduria_entomon_per_mass' (double) with 1,066 unique values and 0% NA
#>         new variable 'prop_saduria2' (double) with 782 unique values and 0% NA

fle_fr2 <- fle_fr2 %>%
  mutate(saduria_entomon_per_mass = saduria_entomon_tot/pred_weight_g,
         prop_saduria2 = ifelse(prop_saduria == 0, prop_saduria + 0.0001, prop_saduria),
         prop_saduria2 = ifelse(prop_saduria == 1, prop_saduria - 0.0001, prop_saduria2)) 
#> mutate: new variable 'saduria_entomon_per_mass' (double) with 430 unique values and 0% NA
#>         new variable 'prop_saduria2' (double) with 300 unique values and 0% NA

# OK, proportion models look absolutely AWFUL. I will model stomach content with a tweedie instead

# Proportion Saduria
# mcodsad <- sdmTMB(
#   data = cod_fr2, 
#   formula = prop_saduria2 ~ 0 + as.factor(year) + predfle_density_sc + predcod_density_sc,
#   time = "year", fields = "IID", spatial_only = FALSE, spde = pcod_spde, family = student(df = 5))

mcodsad <- sdmTMB(
  data = cod_fr2, 
  formula = saduria_entomon_per_mass ~ 0 + as.factor(year) + s(predfle_density_sc) + s(predcod_density_sc),
  time = "year", fields = "IID", spatial_only = FALSE, spde = pcod_spde, family = tweedie())
#> Warning: Detected a `s()` smoother. Smoothers are penalized in sdmTMB as
#> of version 0.0.19, but used to be unpenalized.
#> You no longer need to specify `k`, specifing `k` now indicates
#> an upper limit on the basis functions, and the degree of wiggliness
#> is determined by the data.

cod_fr2$resids_mcodsad <- residuals(mcodsad) # randomized quantile residuals
#> Warning: Detected a `s()` smoother. Smoothers are penalized in sdmTMB as
#> of version 0.0.19, but used to be unpenalized.
#> You no longer need to specify `k`, specifing `k` now indicates
#> an upper limit on the basis functions, and the degree of wiggliness
#> is determined by the data.
hist(cod_fr2$resids_mcodsad)

qqnorm(cod_fr2$resids_mcodsad); abline(a = 0, b = 1)


tidy(mcodsad)
#>                       term   estimate std.error
#> 1      as.factor(year)2007  -7.559057 0.7891209
#> 2      as.factor(year)2008  -7.800466 0.8602433
#> 3      as.factor(year)2009  -8.748290 0.8173057
#> 4      as.factor(year)2010  -9.210658 1.0224320
#> 5      as.factor(year)2011 -10.058870 0.9114980
#> 6      as.factor(year)2012  -8.706427 0.8372932
#> 7      as.factor(year)2013  -9.184026 0.7839526
#> 8      as.factor(year)2015  -8.567464 0.7569644
#> 9      as.factor(year)2016  -9.869242 0.6801726
#> 10     as.factor(year)2017 -10.710367 0.7368018
#> 11     as.factor(year)2018  -9.698110 0.7087329
#> 12     as.factor(year)2019 -10.655589 0.8360574
#> 13 s(predfle_density_sc).1         NA        NA
#> 14 s(predfle_density_sc).2         NA        NA
#> 15 s(predfle_density_sc).3         NA        NA
#> 16 s(predfle_density_sc).4         NA        NA
#> 17 s(predfle_density_sc).5         NA        NA
#> 18 s(predfle_density_sc).6         NA        NA
#> 19 s(predfle_density_sc).7         NA        NA
#> 20 s(predfle_density_sc).8         NA        NA
#> 21 s(predfle_density_sc).9         NA        NA
#> 22 s(predcod_density_sc).1         NA        NA
#> 23 s(predcod_density_sc).2         NA        NA
#> 24 s(predcod_density_sc).3         NA        NA
#> 25 s(predcod_density_sc).4         NA        NA
#> 26 s(predcod_density_sc).5         NA        NA
#> 27 s(predcod_density_sc).6         NA        NA
#> 28 s(predcod_density_sc).7         NA        NA
#> 29 s(predcod_density_sc).8         NA        NA
#> 30 s(predcod_density_sc).9         NA        NA
summary(mcodsad)
#> Spatial model fit by ML ['sdmTMB']
#> Formula: saduria_entomon_per_mass ~ 0 + as.factor(year) + s(predfle_density_sc) + 
#> Formula:     s(predcod_density_sc)
#> Time column: "year"
#> SPDE: pcod_spde
#> Data: cod_fr2
#> Family: tweedie(link = 'log')
#>                     coef.est coef.se
#> as.factor(year)2007    -7.56    0.79
#> as.factor(year)2008    -7.80    0.86
#> as.factor(year)2009    -8.75    0.82
#> as.factor(year)2010    -9.21    1.02
#> as.factor(year)2011   -10.06    0.91
#> as.factor(year)2012    -8.71    0.84
#> as.factor(year)2013    -9.18    0.78
#> as.factor(year)2015    -8.57    0.76
#> as.factor(year)2016    -9.87    0.68
#> as.factor(year)2017   -10.71    0.74
#> as.factor(year)2018    -9.70    0.71
#> as.factor(year)2019   -10.66    0.84
#> 
#> Matern range parameter: 42.38
#> Dispersion parameter: 0.17
#> Spatial SD: 2.29
#> Spatiotemporal SD: 1.43
#> ML criterion at convergence: -3251.017
#> 
#> See ?tidy.sdmTMB to extract these values as a data frame.

mflesad <- sdmTMB(
  data = fle_fr2, 
  formula = saduria_entomon_per_mass ~ 0 + as.factor(year) + s(predfle_density_sc) + s(predcod_density_sc),
  time = "year", fields = "IID", spatial_only = FALSE, spde = pfle_spde, family = tweedie())
#> Warning: Detected a `s()` smoother. Smoothers are penalized in sdmTMB as
#> of version 0.0.19, but used to be unpenalized.
#> You no longer need to specify `k`, specifing `k` now indicates
#> an upper limit on the basis functions, and the degree of wiggliness
#> is determined by the data.

tidy(mflesad)
#>                       term  estimate std.error
#> 1      as.factor(year)2015 -7.324439  1.072835
#> 2      as.factor(year)2016 -7.231370  1.020008
#> 3      as.factor(year)2017 -7.521263  1.017771
#> 4      as.factor(year)2018 -7.248929  1.060067
#> 5      as.factor(year)2019 -7.649248  1.130401
#> 6  s(predfle_density_sc).1        NA        NA
#> 7  s(predfle_density_sc).2        NA        NA
#> 8  s(predfle_density_sc).3        NA        NA
#> 9  s(predfle_density_sc).4        NA        NA
#> 10 s(predfle_density_sc).5        NA        NA
#> 11 s(predfle_density_sc).6        NA        NA
#> 12 s(predfle_density_sc).7        NA        NA
#> 13 s(predfle_density_sc).8        NA        NA
#> 14 s(predfle_density_sc).9        NA        NA
#> 15 s(predcod_density_sc).1        NA        NA
#> 16 s(predcod_density_sc).2        NA        NA
#> 17 s(predcod_density_sc).3        NA        NA
#> 18 s(predcod_density_sc).4        NA        NA
#> 19 s(predcod_density_sc).5        NA        NA
#> 20 s(predcod_density_sc).6        NA        NA
#> 21 s(predcod_density_sc).7        NA        NA
#> 22 s(predcod_density_sc).8        NA        NA
#> 23 s(predcod_density_sc).9        NA        NA
summary(mflesad)
#> Spatial model fit by ML ['sdmTMB']
#> Formula: saduria_entomon_per_mass ~ 0 + as.factor(year) + s(predfle_density_sc) + 
#> Formula:     s(predcod_density_sc)
#> Time column: "year"
#> SPDE: pfle_spde
#> Data: fle_fr2
#> Family: tweedie(link = 'log')
#>                     coef.est coef.se
#> as.factor(year)2015    -7.32    1.07
#> as.factor(year)2016    -7.23    1.02
#> as.factor(year)2017    -7.52    1.02
#> as.factor(year)2018    -7.25    1.06
#> as.factor(year)2019    -7.65    1.13
#> 
#> Matern range parameter: 79.85
#> Dispersion parameter: 0.22
#> Spatial SD: 2.73
#> Spatiotemporal SD: 0.68
#> ML criterion at convergence: -1155.956
#> 
#> See ?tidy.sdmTMB to extract these values as a data frame.

fle_fr2$resids_mflesad <- residuals(mflesad) # randomized quantile residuals
#> Warning: Detected a `s()` smoother. Smoothers are penalized in sdmTMB as
#> of version 0.0.19, but used to be unpenalized.
#> You no longer need to specify `k`, specifing `k` now indicates
#> an upper limit on the basis functions, and the degree of wiggliness
#> is determined by the data.
hist(fle_fr2$resids_mflesad)

qqnorm(fle_fr2$resids_mflesad); abline(a = 0, b = 1)


# Plot marginal effects
## Flounder
# Flounder density
nd_fle_fle <- data.frame(predfle_density_sc = seq(min(fle_fr2$predfle_density_sc),
                                                  max(fle_fr2$predfle_density_sc), length.out = 100))

nd_fle_fle$year <- 2018L
nd_fle_fle$predcod_density_sc <- 0

pred_fle_fle <- predict(mflesad, newdata = nd_fle_fle, se_fit = TRUE, re_form = NA)
#> Warning: Detected a `s()` smoother. Smoothers are penalized in sdmTMB as
#> of version 0.0.19, but used to be unpenalized.
#> You no longer need to specify `k`, specifing `k` now indicates
#> an upper limit on the basis functions, and the degree of wiggliness
#> is determined by the data.

p1 <- ggplot(pred_fle_fle, aes(predfle_density_sc, exp(est),
  ymin = exp(est - 1.96 * est_se), ymax = exp(est + 1.96 * est_se))) +
  geom_line() + geom_ribbon(alpha = 0.4) +
  theme_classic(base_size = 14) + 
  ggtitle("Fle stomach, flounder margin")

# Cod density
nd_fle_cod <- data.frame(predcod_density_sc = seq(min(fle_fr2$predcod_density_sc),
                                                  max(fle_fr2$predcod_density_sc), length.out = 100))

nd_fle_cod$year <- 2018L
nd_fle_cod$predfle_density_sc <- 0

pred_fle_cod <- predict(mflesad, newdata = nd_fle_cod, se_fit = TRUE, re_form = NA)
#> Warning: Detected a `s()` smoother. Smoothers are penalized in sdmTMB as
#> of version 0.0.19, but used to be unpenalized.
#> You no longer need to specify `k`, specifing `k` now indicates
#> an upper limit on the basis functions, and the degree of wiggliness
#> is determined by the data.

p2 <- ggplot(pred_fle_cod, aes(predcod_density_sc, exp(est),
  ymin = exp(est - 1.96 * est_se), ymax = exp(est + 1.96 * est_se))) +
  geom_line() + geom_ribbon(alpha = 0.4) + 
  theme_classic(base_size = 14) +
  ggtitle("Fle stomach, cod margin")

flemargin <- p1 + p2 + plot_annotation(title = "Flounder marginal effects",
                                       theme = theme(plot.title = element_text(size = 16)))

## Cod
# Flounder density
nd_cod_fle <- data.frame(predfle_density_sc = seq(min(cod_fr2$predfle_density_sc),
                                                  max(cod_fr2$predfle_density_sc), length.out = 100))

nd_cod_fle$year <- 2018L
nd_cod_fle$predcod_density_sc <- 0

pred_cod_fle <- predict(mcodsad, newdata = nd_cod_fle, se_fit = TRUE, re_form = NA)
#> Warning: Detected a `s()` smoother. Smoothers are penalized in sdmTMB as
#> of version 0.0.19, but used to be unpenalized.
#> You no longer need to specify `k`, specifing `k` now indicates
#> an upper limit on the basis functions, and the degree of wiggliness
#> is determined by the data.

p3 <- ggplot(pred_cod_fle, aes(predfle_density_sc, exp(est),
  ymin = exp(est - 1.96 * est_se), ymax = exp(est + 1.96 * est_se))) +
  geom_line() + geom_ribbon(alpha = 0.4) +
  theme_classic(base_size = 14) +
  ggtitle("Cod stomach, flounder margin")

# Cod density
nd_cod_cod <- data.frame(predcod_density_sc = seq(min(cod_fr2$predcod_density_sc),
                                                  max(cod_fr2$predcod_density_sc), length.out = 100))

nd_cod_cod$year <- 2018L
nd_cod_cod$predfle_density_sc <- 0

pred_cod_cod <- predict(mcodsad, newdata = nd_cod_cod, se_fit = TRUE, re_form = NA)
#> Warning: Detected a `s()` smoother. Smoothers are penalized in sdmTMB as
#> of version 0.0.19, but used to be unpenalized.
#> You no longer need to specify `k`, specifing `k` now indicates
#> an upper limit on the basis functions, and the degree of wiggliness
#> is determined by the data.

p4 <- ggplot(pred_cod_cod, aes(predcod_density_sc, exp(est),
  ymin = exp(est - 1.96 * est_se), ymax = exp(est + 1.96 * est_se))) +
  geom_line() + geom_ribbon(alpha = 0.4) + 
  theme_classic(base_size = 14) +
  ggtitle("Cod stomach, cod margin")


(p1 | p2) / (p3 | p4 )